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Summary 

A detailed analysis ot the accuracy of several techniques 
recently developed for integrating stiff ordinary differential 
equations is presented. The techniques include two general- 
purpose codes, EPISODE and LSODE, which were developed 
for an arbitrary system of ordinary differential equations, and 
three specialized codes, CHEMEQ, CREK1D, and GCKP84, 
which were developed specifically to solve chemical kinetic 
rate equations. The accuracy study is made by applying these 
codes to two practical combustion kinetics problems. Each 
problem describes adiabatic, homogeneous, gas-phase chemical 
reactions at constant pressure and includes all three combustion 
regimes: induction, heat release, and equilibration. To 
illustrate the error variation in the different combustion 
regimes, the species are divided into three types: reactants, 
intermediates, and products. Error versus time plots are 
presented for each species type and the temperature. These 
plots show that CHEMEQ is the most accurate code during 
induction and early heat release. During late heat release and 
equilibration, however, the other codes are more accurate. A 
single global quantity, a mean integrated root-mean- square 
error, that measures the average error incurred in solving the 
complete problem is used to compare the accuracy of the 
codes. Among the codes examined, LSODE is the most 
accurate for solving chemical kinetics problems. It is also the 
most efficient code, in the sense that it requires the least 
computational work to attain a specified accuracy. An 
important finding is that use of the algebraic enthalpy 
conservation equation to compute the temperature can be 
more accurate and efficient than integrating the temperature 
differential equation. 


Introduction 

Many practical problems arising in chemically reacting flows 
require the simultaneous solution of large sets of coupled 
ordinary differential equations (ODEs) which describe the 
time rate of change of chemical species concentrations and 
temperature. Examples of such problems include the develop- 
ment and validation of reaction mechanisms, combustion ot 
fuel-air mixtures, and pollutant formation and destruction. 

The main difficulty in using classical methods, such as the 
popular explicit Runge-Kutta method (e.g., ref. 1), to solve 
large sets of chemical kinetic rate equations is that of 


“ stiffness.’ ' The property of stiffness arises in chemical 
kinetics because of the widely varying time constants for 
different species. For free radicals the relaxation time is on 
the order of microseconds, whereas the nitric oxide formation 
time is on the order of milliseconds. To satisfy the stability 
requirements that errors in the numerical solution remain 
bounded as the calculation proceeds in time, classical methods 
must use extremely small step sizes, as illustrated in 
references 2 and 3 for the explicit Runge-Kutta method in 
solving combustion kinetics problems. Consequently, these 
methods require prohibitive amounts of computer time to solve 
a practical chemical kinetics problem. 

Numerous approaches have been proposed for stiff ODE s 
to remove the stability restriction on the step size. In Part I 
of this effort (ref. 2) and other recent publications (refs. 3 to 5), 
several techniques were examined, and detailed comparisons 
of their computational work requirements for solving com- 
bustion kinetic rate equations were made. The methods 
examined in these studies include the general-purpose packages 
EPISODE and LSODE (refs. 6 to 9), which were developed 
for an arbitrary system of ODEs, and the specialized codes 
CHEMEQ (ref. 10), CREK1D (refs. 1 1 to 14), and GCKP84 
(ref. 15), which have all been developed specifically to 
integrate chemical kinetic rate equations. In the present work 
the accuracy of these techniques in solving combustion kinetic 
rate equations is examined. 

In general, numerical methods generate approximate 
solutions to the governing ODE's at discrete points in time. 
To maintain accuracy of the numerical solution, they require 
that the estimated error incurred on each time step be less than 
a user-specified local error tolerance. This result is usually 
achieved by restricting the size of the time step. Some solvers, 
in addition, adjust the order of the numerical approximation 
when appropriate. In either case, only the estimated local error, 
that is, the estimate of the error incurred in advancing the 
numerical solution by one time step, is controlled. However, 
the quantity that is of interest to the user is the global error, 
which is the deviation of the numerical approximation from 
the exact solution and which generally accumulates in a 
nontrivial manner from the local errors. 

In the present paper, a detailed study of the estimated global 
error incurred by the above techniques in solving combustion 
kinetic rate equations is presented. Also presented is a study 
of the variation of the global error with the user-specified local 
tolerance and an examination of the computational cost, 
measured by the required CPU execution time, associated with 
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attaining desired accuracy. The paper concludes with two 
appendixes: Appendix A describes the methods examined in 
this study, and appendix B describes the procedure used to 
solve the algebraie enthalpy conservation equation for the 
temperature. 


Symbols 


pre-exponential constants in forward and reverse 
rate coefficients for reaction j (eqs. (6) and (7)), 
units depend on reaction type 

ATOL, local absolute error tolerance for / th component, 
required by LSODE (eq. (20)). 

ATOLSP local absolute error tolerance used with LSODE 
for all species mole numbers 

Bj exponent-on-ten in pre-exponential constant for 

forward rate coefficient of reaction j, where B j = 
logio^, arbitrary units 

C G local error test constant used in GCKP84 (eq. (21)) 

c pJ constant-pressure molar specific heat of species 

/. J/kmole K 


estimated local truncation error in component 
at /„ 

cumulative difference between converged and pre- 
dicted values of (dYjdt) at used by GCKP84, 
units depend on component i 

BjyE-j activation energy in forward and reverse rate 
coefficients for reaction j (eqs. (6) and (7)), cal/mole 

£ r ms mean integrated root-mcan-square global error 
(eq. (34)) 


EPS 


ERMAX 

EWT, 


e U 


L rms 




e r 


for EPISODE and GCKP84; local relative error 
tolerance for species with initially nonzero mole 
numbers and temperature, and local absolute 
error tolerance for species with initially zero mole 
numbers; for LSODE: local relative error tolerance 
for all components; for CHEMEQ and GREK 1 D: 
local relative convergence criterion for all 
components. 

relative error tolerance for Newton-Raphson iteration 
for temperature 

local error weight for / th component, used by 
LSODE (eqs. (19) and (20)) 

estimated global error in /' th component (eq. (24)) 
estimated global error for / th species of type j 

root-mean-square norm of the estimated global 
errors for all variables (eq. (27)) 

root-mean-square norm of the estimated global 
errors for species of type j (eq. (26)) 

estimated global error in temperature (eq. (25)) 


3 generalized algorithm for Y ft (eq. (14)) 

f time rate of change of / th component, units depend 

on i 


K 

ho 

IERROR 

ITOL 

k r k -j 

N 

W-J 

Hr 


initial mixture mass-specific enthalpy, J/kg 

molar-specific enthalpy of species /, J/kmole 

step size used on the n lh step, s 

initial step length to be attempted by integrator, s 

error control indicator for EPISODE 

error control indicator for LSODE 

forward and reverse rate coefficients for reaction 
j (eqs. (6) and (7)), units depend on reaction type 

total number of first-order ordinary differential 
equations 

temperature exponent in forward and reverse rate 
coefficients for reaction j (eqs. (6) and (7)) 

total number of elementary chemical reactions in 
reaction mechanism 


total number of chemical species in reacting gas 
mixture 


P1,P2 

P 

R 

RTOL 

T 

AT 

T S r 

TINY 

t 

^end 

h) 


test problems I and 2, respectively 
pressure, N/nr 

universal gas constant in cal/mole K 
universal gas constant in J/kmole K 

molar forward and reverse rates per unit volume 
for reaction / (eqs. (4) and (5)), kmole/m* s 
local relative error tolerance for LSODE 
temperature, K 

maximum temperature change allowed before 
reaction rate coefficients and thermodynamic 
properties are updated in CREK1D, K 

standard solution value for temperature, K 

minimum species mole number values allowed in 
CHEMEQ and CREK ID 
reaction time, s 

final time (> 1 ms) at which numerical solution 
is generated, s 

time reached on the /z ,h integration step, s 
initial time, s 


V 

X ( 

*i.st 

'*mm 

Yu 


dm| 


reacting gas velocity, m/s 

chemical symbol for / ,h species 

mole fraction of species i 

standard solution value for mole fraction of species i 

mole fraction value corresponding to <j mjn (eq. (33)) 

numerical solution of the / ,h component at f /n units 
depend on i 

value obtained for Y jn on w lh iteration, units 
depend on i 


2 



y M 
* j'.#i 





V 

1 mux./ 

V/ 


a 


*n 


P 

^inin 


value obtained for K,„ on m* iteration, units 
depend on i 

predicted value of Y it ,„ units depend on / 
numerical solution of the i ,h component at t 
generated by using exact past values, units depend 

on / 

local weight for i' ,h component, used by EPISODE 
and GCKP84 (eqs. (17), (18), and (21)) 
exact solution for the t' h component, units depend 
on i 

constant in generalized algorithm for Y n (eq. (14)) 
global error at t n (eq. (15)) 

stoichiometric coefficients of species i in forward 
and reverse reaction j (eq. (1)K number o 
kilomoles of species i in elementary reaction j as 
a reactant and as a product, respectively 

mixture mass density, kg/m 

mole number of species /. kmole species i/kg mixture 
mole number value at which local error control in 
LSODE is equally relative and absolute (eq. (- *•)) 


Governing Differential and 
Algebraic Equations 

The ordinary differential equations describing homogeneous 
gas-phase chemical reactions ot the type 


N s 

£ D tjX; 7=1 N* 

1=1 i = 1 


( 1 ) 


Ns 


Rj = k < n (pa,) '" 


(4) 


/=1 


and 




(5) 


/=1 


where the forward (k,) and reverse (*_;) rate coefficients are 
given by the modified Arrhenius expressions: 


A, = A,T n • exp (-Ej/RT) 
k_j = A_jT s i exp ( — E-j/RT) 


( 6 ) 

(7) 


/n {1 \ „ -ind v are the stoichiometric 

In equations (1) to (/), Vjj anu v ij 

coefficients of species / (with chemical symbol X,) in reaction 

j as a reactant and as a product, respectively; N 5 i s the tola 

number of distinct chemical species (reacting and inert) 

gas mixture; N R is the total number of independent reactions 

in the mechanism; a, is the mole number of spec.es Mm 

kilomoles of species / per kilogram of mixture); t is the time 

(in seconds); p is the gas mixture mass density (in 

per cubic meter); T is the temperature (in kelvins) R s the 

universal gas constant (in calories per mole per kelvm) and 

A 4 N N-j E r and E .j are constants in the modified 

Arrhenius ^ expressions for *, and The reverse rate 

coefficient parameters arc calculated from the forward 

coefficient parameters and the concentration equilibrium 

constants b, using the principle of detailed balanemg (ref. 16b 

In this paper, as in the companion paper (ref. -). 

is restricted to adiabatic, constant-pressure chemical reactions. 
For ' such problems, the following enthalpy conservation 
equation constitutes an algebraic constraint on equations (2) 


are as follows: 


c !^=f i (a k ,T) i,k = 1 N, 

dt 

o,(f 0 ) = g' ven 
T{t 0 ) = given 


( 2 ) 


where /, the total formation rate of species /. is given by 


a,!), = H 0 , 

i = i 


( 8 ) 


■e ft is the molar-specific enthalpy of species i (in joules 
kilomole) and H 0 is the initial mixture mass-specific 
alpy (in joules per kilogram). Equation (8) can be 
rentiated with respect to time to give the following OD 
he temperature 


7=1 

The molar reaction rates per unit volume, Rj and are 

given by the law of mass action (e.g., ref. 16): 


dT 

dt 


i=i 

N, 

/ = 1 


(9) 
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where c p , is the constant-pressure molar-specific heat of 
species / (in joules per kilomole per kelvin). Either equation 
(8) or (9) can he included in the equation set. We explore the 
use of both these equations and examine their effects on 
solution accuracy and computational cost. 

The mass density of the mixture is given by the ideal gas 
equation of state 


P=p/(RJo m ) (10) 

where/? is the absolute pressure (in newtons per square meter), 
R t , is the universal gas constant (in joules per kilomole per 
kelvin), and the reciprocal of the mean molar mass of the 
mixture, is given by 


N x 

a >» = E °> do 

/=] 

For constant pressure problems the following density ODE 
can be obtained from equation (10) by differentiating it with 
respect to time and then rearranging terms in the resulting 
expression: 7 


dp 

dt 


7 dt r* dt) 


( 12 ) 


Either equation ( 1 0) or ( 1 2) can be used to compute the density . 
The code GCKP84, which allows the pressure to vary, solves 
forp by integrating its ODE (eq. (12)). With the other codes, 
however, we obtain p by using equation (10). Indeed, p is 
implicitly replaced by the right-hand side of equation (10) and 
does not appear as a variable. We, therefore, exclude density 
from our discussion, including statement of the problem, and 
restrict attention to solving for the other N s + 1 quantities. 


The packages EPISODE and LSODE consist of a variable- 
order, variable-step implicit Adams method (suitable for 
nonstiff problems) and a variable-order, variable-step backward 
differentiation formula method (suitable for stiff problems: 
e.g., refs. 1 and 17). Both methods use a standard predictor 
and a variety of corrector formulas— from functional iteration 
to a modified Newton iteration— is included. The Jacobian 
rnatri x df/dy, where y is the vector of dependent variables 
and/ - dy/dt, is computed either numerically or with a user- 
supplied subroutine. In part I of this investigation (ref. 2) all 
options relevant to the problem of chemical kinetics were 
attempted, and the stiff method with Newton iteration and user- 
supplied analytical Jacobian matrix was found to be the fastest. 
Therefore, only this option is used in examining the accuracy 
of EPISODE and LSODE. 

The general chemical kinetics program GCKP84 uses the 
integration technique developed by Zeleznik and McBride 
(ref. 18). The algorithm is essentially a revised version of the 
GEAR package (ref. 19), which contains the same two integration 
methods as EPISODE and LSODE and several iteration 
techniques. For reasons given above we restrict attention to 
the stiff method with Newton iteration using an analytical 
Jacobian matrix. GCKP84 includes corrective actions if the 
physically impossible situation of negative concentrations, 
temperature, density, or velocity arises. 

In CHEMEQ, at the start of each time step the ODE's are 
separated into two classes: stiff and normal. For equations 
classified as normal, a classical second-order predictor-corrector 
method, the trapezoidal rule, is used. For the stiff equations 
a simple stable asymptotic integration formula is used. 

The code CREK1D is based on the exponentially fitted 
trapezoidal rule developed by Liniger and Willoughby (ref 20) 
and Brandon (refs. 21 and 22). This code includes special 
treatment of ill-posed initial conditions and automatic selection 
of Jacobi-Newton iteration or Newton iteration. 


Problem Statement 

The initial value problem may be stated as follows: Given 
(1) at time t = t 0 , values for the species mole numbers, a, 

(/ = 1 N s)< an d the temperature, T, (2) the pressure, p, 

and (3) the reaction mechanism, find, at the end of a prescribed 
time interval, the mixture composition and temperature. 

Methods and Codes Examined 

The codes examined in this study include the general-purpose 
packages EPISODE and LSODE (refs. 6 to 9) and the 
specialized techniques CHEMEQ (ref. 10), CREK1D (refs. 

1 1 to 14), and GCKP84 (ref. 15). The methods used in these 
codes are summarized below and are discussed in detail in 
appendix A. 


Error Considerations 


ui me uuicrciH codes 

examined are discussed. In general, numerical methods replace 
the differential equations with difference equations and solve 
them step by step. Starting with the known initial conditions 
y(/ 0 ) at /„ numerical approximations Y„ to the exact solution 
y(f„) of the ODE’s are generated at discrete points in time 
( { n ( n — L2,...)), until the end of the integration interval is 
reached. At each t„ the numerical method provides a rule for 
generating the approximate solution Y„ in terms of computed 
quantities at one or more previous times. 

For the scalar differential equation. 


dy 

dt 


= /(v) 


y(/ 0 ) = given 


(13) 
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the algorithms used for Y„ in all of the codes can be written as 
Y„ = cxYn - 1 + h’t $ (h n ,Y n ’Jii’Yii-\'fn- ^ 

where a is a constant, r„_, is the approximate solution at f,-,, 
h„ (equal to t„ - f„-,) is the step size used on the step 
[/, „f„l and f„ -j =f(Y„ _,)• Because equation (14) involves 
the unknown quantity Y„ its solution generally requires an 
iterative procedure. Starting with the predicted value (an initial 
guess), denoted by if 1 , improved estimates if (m - 
arc generated until the iteration converges, that is. until the 
difference in two successive approximations approaches zero 
within a specified accuracy. 

During the calculation procedure, errors called discretization 
or truncation errors are introduced into the numerical solution 
because of the approximation of the ODEN by difference 
equations. Two measures can be defined for this error, w hi c 
is a property of the numerical method (e.g.. ref. 2d). e 
alobal discretization error e„ at any t„ is the difference 
between the computed approximation Y n and the exact 
solution y (/„).* 

e„ = Y„ - y(t„) (l5) 

It is the quantity that the user wants to know and control. The 
local truncation error d„ at t„ is the error in the numerical 
approximation Y„ that is generated on the step Un-i> f „l y 
using exact past values: 


It is the quantity that ODE solvers generally control. The two 
discretization errors are illustrated in figure 1 for a sing e 
ODE. 

The codes examined in this study require the user to specify 
values for one or more local tolerance parameters, which 
control the accuracy of the numerical solution. Now as 
discussed below, the same error control is not used in all codes. 
Nevertheless, for convenience, for all codes the same notation, 
EPS is used to denote the local tolerance quantity, or the 
primary one if several are required. In EPISODE the local 
truncation error vector d„ satisfies the inequality 


< EPS 


where N is the number of ODE’s, is the estimated local 
truncation error in the i' ,h component at t„, and for the error 
control used 


Exact solution 

#,■ Numerical solution 


f.nux, = maX W Y i,n~\\ ’ 1^1-2 1 1 

= max [1 , 


n rf 

1‘iC 


y(t 0 ), 


Figure 1 . — Numerical solutions and truncations error t>pcs for the sinj.le 
dxldt = f(y) The exact solution to the ODE is denoted by v(r) t 
numerical solutions obtained with the initial condition = > (»o) ■« 
denoted by solid circles. The solid square denotes the numerical solution 
obtained at t„ by using exact past values. The iocal .runcation crror is 
denoted by d. the global truncation error by t. and the step length b. . 

where the vertical bars denote absolute value. The error control 
selected to be performed by LSODE is given by 


where 


EWT, = EPS [y,., I -il + ATOL, 


where ATOL, is the user-supplied local absolute error 

tolerance for the i ,h component. 

In GCKP84 the local error test satisfies the inequality 


; 5 ic 


< Cc. EPS 


for i that satisfy y,(fo) ^ 0, 

for / that satisfy y,(fo) = 0, 
(18) 


where, E, „ contains the cumulative difference between t e 
converged and predicted values of the derivative (dYJdt) X. 
t and where C G is a constant. The quantity T m ax./ has t e 
same meaning as in EPISODE (see eq. 18)). 

The codes CHEMEQ and CREK1D do not control the 
estimated local truncation error. The solution is accepted w en 
the magnitude of the normalized difference in successive 
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eslimaics S'") is less ,ha„ a specified amomu. 

rherefore, these codes control only the error in the solution 
to the difference equation of the method. In CHEMEQ each 
component i satisfies the inequality 

j V.lw + 1 1 yjwj 1 

1 f./l I 

min finir +l| j,|nl;'"l! ~ EPS (22) 

The convergence criterion used in CREK 1 D is given by 

(1 y (Y}:r n -Yw v\ 2 

y-Vy ™ \ yjml J J — EPS (23) 

It is clear from the above discussion that the user-supplied 
ocal tolerance EPS does not have the same meaning for all 
codes. In LSODE it is the local relative error tolerance for 
all variables and is a measure of the number of accurate 
significant figures in the numerical solution. In EPISODE and 
GCKP84, however, as discussed in the section “Computational 
Procedure, EPS is the local relative error tolerance for only 
variables with nonzero initial values, such as the temperature. 
For species with zero initial mole numbers EPS is the local 
absolute error tolerance and is a measure of the largest number 
that may be neglected. In contrast to these three codes. 
CHEMEQ and CREK ID do not control the local truncation 
error, and EPS is the local relative convergence criterion, or 
error m the solution to the difference equation. However as 
described in appendix A, although CREK 1 D does not test that 
the estimated local truncation error is within a prescribed 
bound, the step length calculation procedure attempts to 
achieve this result. The step length to be attempted next is 
selected such that the current estimate of the local truncation 
error normalized by the solution is at most equal to EPS. 
Because of these differences in the meanings of EPS it will 
be referred to as simply the local tolerance. 

Evaluation of Temperature 

Of the codes examined in the present study, only GCKP84 
and CREK ID were written explicitly for nonisothermal chemical 
reactions. These methods, therefore, have built-in procedures 
for calculating the temperature. For the other codes however 
the temperature has to be calculated along with the mixture 
composition. In the present study (as in ref. 2), the temperature 
was computed using one of the two methods outlined below 
In method A the temperature was calculated from the initial 
mixture mass-specific enthalpy H 0 and the solution for the 
species mole numbers returned by the integrator by using the 
algebraic enthalpy conservation equation (8). This equation 
was solved for the temperature by using a Newton-Raphson 
iterative technique, with a user-supplied local relative error 
tolerance, ERMAX (as described in appendix B). In this 
method, the temperature is not an explicit dependent variable. 


so the number of ODE's is equal to the number ( N s ) of 
species and the Jacobian matrix is of size A; x N s The 
integrator, therefore, tracks only the solution for the species 
mole numbers. The temperature was also computed when the 
species time derivatives and the Jacobian matrix were evaluated. 

In method B the temperature was treated as an additional 
dependent variable and evaluated by solving its ODE (eq (9)) 
In this method, the number of ODE’s is equal to N s + I , the 
Jacobian matrix is of size (N s + 1) x (N s + 1), and the inte- 
grator tracks the solutions for both the species mole numbers 
and the temperature. 

The following naming convention was adopted. Techniques 
using method A were given the suffix A (EPISODE- A. etc.), 
and those using method B were given the suffix B (EPISODE- 
B, etc.). 


uaiiMci oeiween tne 

reacting gas mixture and its surroundings and must therefore 
use an ODE to solve for the temperature. It also includes the 
density and velocity, V. of the gas mixture as dependent variables 
and evaluates them by integrating their ODE’s. (For the static 
test problems used in this study the velocity ODE is given 
trivially by dV/dt = 0, V(t 0 ) = 0.) Consequently, the number 
of ODE s solved by GCKP84 is equal to N s + 3, and the 
Jacobian matrix is of size (N s + 3) x (/V9 + 3) 

CREK ID computes the temperature by solving the algebraic- 
enthalpy conservation equation (8). However, the calculation 
procedure is different from that used in method A. In CREK I D 
the mixed differential-algebraic system of equations (2) and 
( ) is solved simultaneously, whereas method A solves 
equation (8) after the species ODE’s have been integrated 
over a time step. Thus, although the number of ODE’s solved 

by CREK ID is equal to N s , the Jacobian matrix is of size 
(Vs + 1) x (N s + ]). 


Test Problems 

The algorithms examined in the present study were applied 
to the same two test problems used in our previous work 
(ref. 2). Both problems describe adiabatic, constant pressure 
transient, batch chemical reactions and include all three com- 
bustion regimes: induction, heat release, and equilibration. 

Test problem 1 describes the ignition and subsequent combus- 
tion of a mixture of 33 percent carbon monoxide and 67 percent 
hydrogen with 100 percent theoretical air at an initial temperature 
of 1000 K and a pressure of 10 atm. It comprises 12 reactions 
which describe the temporal evolution of 1 1 reacting species 
(CO, C0 2 , H, H 2 , H 2 0, N, NO, N 2 , o, OH. and 0 2 ). Test 
problem 2 describes the ignition and subsequent combustion 
of a stoichiometric hydrogen-air mixture at a pressure of 2 atm 
and an initial temperature of 1500 K. It involves 30 reactions 
among 15 species (Ar, C0 2 , H, HO-,, H,, H.O H^O, N 
NO, N0 2 , N 2 , N 2 0, O, OH, and 0 2 ), of which two (Arand 
C0 2 ) are inert. The reaction mechanisms and forward rate 
coefficient parameters for the two test problems are given in 
tables I and II. 
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TABLE I. -REACTION MECHANISM AND FORWARD RATE 
COEFFICIENT PARAMETERS USED FOR TEST PROBLEM 1 
| Rate coefficient k f = 10*' 7*' exp< -E } /RT) .\ 


Reaction 

Reaction 

Rate coefficient parameters 

number. 





./ 



»i 

E , 

keal/mole 

1 

CO + OH = CO, + H 

11.49 

0 

0.596 

2 

H + O, = O + OH 

14.34 


16.492 

3 

H, + 6 = H + OH 

13.48 


9.339 

4 

H,0 + O = OH + OH 

13.92 


18.121 

5 

H + H,0 = H, + OH 

14.0 


19.870 

6 

N + 0 2 = NO + O 

9.81 

1.0 

6.250 

7 

N 2 + o = N + NO 

13.85 

0 

75.506 

8 

NO + M = N + 0 + M 

20.60 

-1.5 

149.025 

9 

H + H + M = H, + M 

18.00 

-1.0 

0 

10 

o + o + m = o 2 + m 

18.14 

-1.0 

0.340 

1 l 

H + OH + M = H,0 + M 

23.88 

-2.6 

0 

12 

H, + O, = OH + OH 

13.00 

0 

43.000 


TABLE IL — REACTION MECHANISM AND FORWARD RATE 
COEFFICIENT PARAMETERS USED FOR TEST PROBLEM 2 
(Rate coefficient k f = 10 fl '7 A 'exp {-Ej/RT).] 


Reaction 

number, 

j 


Rate coefficient parameters 


H + O, = OH + O 14.342 

O + H, = OH + H 10.255 

H : + OH = H,0 + H 13.716 

OH + OH = 6 + H 2 0 12.799 

H + 0 2 + M = HO, + M 15.176 

0 + 0 + M=0, + M 13.756 

H + H + M = H 2 + M 17.919 

H + OH + M = H,0 + M 21.924 

H, + HO, = H : 0 + OH 1 1.857 

H 2 0, + M = OH 4- OH + M 17.068 

H 2 +" 0 2 = OH + OH 13.000 

H + HO, = OH + OH 14.398 

O + HO : - OH + 0 2 13.699 

OH + HO, = H 2 0 Hh 0 2 13.699 

HO, + HO : = H 2 0 2 + 0 2 12.255 

OH + H 2 0 2 = H 2 0 + HO, 13.000 
o + H 2 0 2 = OH + HO, 13.903 

H + H : 0 2 = H : 0 + OH 14.505 

HO, + NO = NO, + OH 13.079 

O + NO, = NO + O, 13.000 

NO T O + M = NO, + M 15.750 
NO, + H = NO + OH 14.462 

N + O, - NO + O 9.806 

O + N 2 = NO + N 14.255 

N + OH = NO + H 13.602 

N,0 + M = N, + O + M 14.152 

O + N,0 = N, + O, 13.794 

O + jsj 2 0 = NO 4- NO 13.491 

N -F NO, = NO + NO 12.556 

OH + N, = N,0 + H 12.505 





Figure 2.— Variation with reaction time of temperature and species mole 
fractions for test problem 1 . 


Figures 2 and 3 present the variations with time of the 
chemical species mole fractions and temperature for test 
problems 1 and 2, respectively. These solutions were generated 
with LSODE-B using a small value (10 ~ 5 ) for the local 
relative error tolerance. Both test problems were integrated 
over a time interval of 1 ms in order to obtain near- 
equilibration of all chemical species and the temperature. 


Computational Procedure 

For each method, global errors in solutions generated with 
a certain value for the local tolerance EPS were estimated by 
comparing them with results obtained with the same method 
and a reduced tolerance. The solutions used as a basis for 
comparison were the most accurate generated and are referred 
to as standard solutions. For example, for CREK1D solutions 
used as standards were generated with CREK1D and 
EPS - 10 \ These standard solutions were used to estimate 
the global errors in results produced with CREK1D and 



Reaction time, t, s 

Figure 3. — Variation with reaction time of temperature and species mole 
fractions for test problem 2. 


EPS = 10 \ 10 3 and 10 4 . The above procedure for 
estimating global errors is reliable provided the technique is 
effective in the sense that reducing the local tolerance actually 
reduces the global error (ref. 24). In any case, in the absence 
of exact solutions the only method for assessing the accuracy 
of an algorithm is to compare the solutions that it produces 
with those obtained with a reduced local tolerance using either 
the same algorithm or a different one. The use of solutions 
generated by each technique as a standard of comparison only 
for itself ensures that the accuracy comparison is not biased 
in favor of any one method or code. 

A typical computational run was performed by first initializing 
the time (f, set equal to zero), species mole numbers, and 
temperature. The integrator was then called with values for 
the necessary input parameters, including the local tolerance 
and the elapsed time (equal to 1 ms for both problems) at which 
the integration was to be terminated. After each step successfully 
executed by the integrator, the current time and the solutions 
for the species mole fractions and temperature were saved. 
This procedure was repeated until the time reached by the 
integrator was greater than or equal to 1 ms. The saved time 
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Temperature, T, K 





values served as input data for the output stations at which 
the standard solution was to be generated. At each ot these 
discrete time values global errors in the species mole fractions 
and temperature were estimated by comparisons with t e 
standard solutions as follows: 


e,U) = 


Xj.Sli 1 ) 


e-fit) 


i = l N, 


TO) 

TsrO) 


(24) 

(25) 


where e,(t) and e r 0) are. respectively, the estimated globa 
errors in the mole fraction jr,(f) of species i and the tem- 
perature T(t) at time t and where x lS rO) and T st 0 ) arc, 
respectively, the standard solution values for the mole fraction 
of species / and the temperature at time t. To prevent the 
possibility of requiring accuracy in species with immeasurably 
small concentrations, global errors were not measure o 
species whose standard solution mole fractions were less than 
0 1 ppm. For such species e,0) was set equal to zero. In this 
way time histories of the global errors in species mole fractions 

and temperature were generated. 

For each technique, standard solutions were generate wi 
a small value for EPS. In addition to EPS and the elapsed time 
at which the integration was to be terminated, other input 
parameters were required by all codes examined. In this paper 
only those input parameters that affect the accuracy o cac 
code are discussed. A more detailed discussion of these 
parameters can be found in part I (ref. 2). 

The user-supplied parameters relevant to solution accuracy 
that are required by LSODE are the error control flag 1TOL, 
which indicates the type of local error control to be performed, 
and the local relative, RTOL, and absolute, ATOL error 
tolerances. Both RTOL and ATOL can be specified either as 
(1) a scalar, so that the same local error tolerance is used tor 
all variables, or (2) an array, so that different values of the 
local error tolerance are used for different variables. In the 
present work the error control given by 1TOL - 2 (for scalar 
RTOL (equal to EPS) and array ATOL, see appendix A) was 
used for reasons given below. Since the same number of 
accurate significant figures is acceptable for all solution 
components, RTOL was specified as a scalar. Now, forth 
test problems examined in this study, the species mole fractions 
and temperature vary widely (figs. 2 and 3), so relative error 
control is appropriate and is the reason for designating the local 
relative error tolerance as the primary tolerance (eq. (20)). 
Pure relative error control can be achieved by specifying a 
value of zero for the local absolute error tolerances. However, 
since many of the species had zero initial concentrations, pure 
relative error control could not be used. To make the error 
control mostly relative, small values were specified for the 
absolute error tolerances for the species mole numbers, or 
convenience the same value (equal to ATOLSP) was used tor 


all species. Since the temperature can never be zero, pure 
relative error control was used for this variable, that is, t e 
local absolute error tolerance for temperature was set equal 
to zero Thus, ATOL was specified as an array. 

The values used for ATOLSP were those obtained in part 
I of ihis smdy (ref. 2) for LSODE-B a s follows: W»h 
EPS = 10 " 5 ATOLSP was progressively decreased until the 
temperature-time trace showed essentially no change with a 
further decrease. The values obtained for ATOLSP by using 
this procedure were 10 14 and 10 , rcspetmc j 4 

problems 1 and 2. For consistency, the same EPS and ATOLSP 
were used with LSODE- A. They were, however, checked to 
ensure that reductions in ATOLSP resulted in essentially the 
same solutions. For reasons given in the next section ERMA 

was set equal to EPS. . 

To make accuracy comparisons among the codes meaningful, 

the same value of EPS (i.e., 10 s ) was used to ^a'c 
standard solutions for GCKP84, CHEMEQ-A, CHEMEQ- 
B. and CREK1D. With EPISODE, however, an EPS value 
of l(T 6 was used because larger values produced physical y 
meaningless results for test problem 1 -little or no change from 
initial values after an elapsed time ot 1 ms. The error contro 
to be performed by this code is selected by means of the Hag 
IERROR (appendix A). For the reasons given above, pure 
relative error control (option IERROR = 2) could not be used 
and the option IERROR = 3 was used, instead. This error 
control is semirelative (see eqs. (17) and (18)). It is relative 
for a variable that is initially nonzero. But for a variable that 
is initially zero, it is absolute until the variable reaches unity 
in magnitude, when it becomes relative. Since none of the mole 
numbers attains a value of unity, the error control is always 
absolute for species with zero initial mole numbers. 

The solution generated with EPISODE depended on the 
value specified for the initial step length (h Q ) to be attempted 
by the integrator. In generating standard so'ut'ons withthis 
code, /i 0 was progressively decreased (with EPS - 10 ) 

until the temperature-time trace showed essentially no change 
with a further decrease. The values obtained for h n by using 
this procedure were 10“’ and 10 * s, res pective y . for test 
problems 1 and 2, for both EPISODE-A and EPISODE-B. 
However, an h 0 value of 10 " 9 s was used for test problem 
2 because it resulted in smaller execution times, as shown in 
table III. For EPISODE-A the savings were modest, but tor 
EPISODE-B they were significant. 

GCKP84 uses the same error control as that selected to e 
performed by EPISODE. It also requires the user to specify 
h n . Since details of the integration technique used in GCKPS4 
were not known, a default value of /i 0 - 10 s wasuse in 
our previous work (refs. 2 to 5). However, Bit.ker and Scull.n 
(ref. 15) have since then set the default value for h 0 at 
5x 10 -8 s. Nevertheless an /i () value of 10 s was used in 
this study to be consistent with part 1 (ref. 2). In addition, as 
shown in the next section, the 10 ft value generally produced 
more accurate results than the new default value, while requiring 
comparable execution times for all EPS used in this study. 
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TABLE III.— EFFECTS OF INITIAL STEP 
LENGTH ON EXECUTION TIMES REQUIRED 
BY EPISODE-A AND B (EPS = 10 h ) 
FOR TEST PROBLEM 2 


Method 

Initial step 
length, 
h 0 , 

s 

CPU 

execution 

time, 

s 

EPISODE A 

10 « 

3.1 


10 g 

3.0 

EPISODE B 

10-* 

14 


10-9 

7.8 J 


In contrast to EPISODE and GCKP84, the other codes 
automatically compute the h {) value to be attempted by the 
integrator. In LSODE the calculation procedure for h {) employs 
the user-specified values for the first output station and the 
local error tolerances. The computed initial step length can 
have an adverse effect on both the computational work and 
the solution generated by the code (ref. 3). The calculation 
procedures used for h 0 in CHEMEQ and CREK I D are based 
on the problem physics (see appendix A) and the computed 
Ao did not cause the above difficulties. 

Both CHEMEQ and CREK ID use a relative convergence 
criterion (eqs. (22) and (23)). The difficulty of applying the 
test when the solution vanishes is avoided by setting mole 
numbers less than a suitably small value, TINY, to be equal 
to TINY. In this study a value for TINY of 10~ 20 was used. 
The only user-specified parameter required by CREK ID that 
affects its accuracy is AT, which is the maximum temperature 
change allowed before the reaction rate coefficients and the 
thermodynamic properties ft, and c p j are updated. Use of this 
parameter increases the efficiency of numerical techniques in 
solving combustion kinetic rate equations (refs. 2 and 3). To 
ensure that the most accurate solutions were used as standards, 
a value of AT = 0 K was used for both test problems. 


Results and Discussion 


make accuracy comparisons between the two temperature 
calculation methods meaningful, ERMAX was set equal to 
EPS, thereby imposing the same local accuracy requirements 
on both methods. Thus, both methods A and B used the same 
error control (i.e., pure relative) and the same local tolerance. 

To facilitate accuracy comparisons among the different 
techniques, the species were divided into three types: reactants 
(R). intermediates (/), and products (P). At each discrete time 
at which global errors had been computed, root-mean-square 
(rms) errors. e rnK j(r) (j = R.I.P ), were computed for all 
three species types as follows: 

/ A* \l /2 

Gms,/ (0 - ( — efj(t) j j = R, /, p (26) 

where e,j(t) ( i = 1 Nj) is the estimated global error at 

time t in the mole fraction of the /' h species of type / and 
where Nj is the number of species of type j. The values of 
Nj (and the species that comprise each subset) are as follows- 
For test problem 1, N H = 4 (CO, H,, N,, and O,), N, = 4 
(H, N, O, and OH), and N P = 3 (CO,. H,0 and NO) For 
test problem 2, N R = 5 (Ar, CO,, H,, N,. and O,), N, = 6 
(H, H0 2 , H,0,, N, O, and OH), and N P = 4 (H,0, NO, 
NO,, and N,0). Although Ar and CO, are inert species, so 
that their mole numbers do not change during the course of 
the reaction, they are classified as reactant species because 
they participate in three-body reactions as catalysts and their 
concentrations affect the rates of these reactions. 

In addition to the rms error for each species type, a single 
rms error for all variables, e rn Jt), was computed at each time 
t by using 




ej(t) + e 2 r (t) 


Gms ( 0 — 


/ = l 


N s + I 


(27) 


The procedure described in the previous section was used 
to study the global errors incurred by the different techniques 
in solving the two test problems. All results presented herein 
were generated on the NASA Lewis Research Center's IBM 
370/3033 computer using single-precision accuracy, except 
GCKP84 which uses double-precision accuracy. 

Both temperature calculation methods A and B were 
attempted with EPISODE, LSODE, and CHEMEQ. The error 
control used in method A is pure relative, and the local relative 
error tolerance is equal to ERMAX (see appendix B) In 
EPISODE-B and CHEMEQ-B the error or convergence 
control for the temperature is pure relative and the local relative 
tolerance is equal to EPS (eqs. (17), (18), and (22)). For 
reasons given previously in the section "Computational 
Procedure" the above remarks apply to LSODE-B also. To 


Figures 4 to 9 present the variations with time of the percent 
rms error in reactants, intermediates, and products and the 
percent error in temperature for test problem I . Similar infor- 
mation is presented for test problem 2 in figures 10 to 17. For 
brevity, test problems 1 and 2 are hereinafter referred to as 
PI and P2, respectively. Note that for clarity the actual errors 
have been magnified in some of the figures. The maximum 
percent rms errors incurred and the reaction times at which 
they occurred are given in tables IV and VI, along with the 
values used for the input parameters discussed in the previous 
section. For each code (except GCKP84) and EPS, these input 
parameters, obtained in part I of this study (ref. 2) by a trial- 
and-error procedure, minimized the execution time required 
to solve the problem. To prevent the possibility of generating 
physically meaningless results by using too large a value of 
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Figure 5.— Variation with time of percent 
I results generated using LSODE-B. 
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Figure 9. -Variation with time of percent mis error in (a) reactants, (b) intermediates, and (c) products, and of(d) percent error in temperature. Test problem 
I results generated using CREK1D. 
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TABLE IV. -ROOT-MEAN-SQUARE ERRORS FOR TEST PROBLEM 1 
(Numbers in parentheses designate powers of 10 .) 
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TABLE VII. -ERRORS FOR TEST PROBLEM 2 


Error in 
temperature 

Standard 

solution. 

K 

2908.5 
1569 9 
1572 4 

2457.7 

1562.1 

1570.3 

1583.7 

1585.9 

1584.4 

1571.1 

2170.9 

1581.1 

2679.5 

1554.8 

1562.2 

1584.8 

2829.8 

2888.6 

2893.4 

2905.5 

2903.2 

2903.2 

2760.4 

2811 .8 

2870.2 
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ATOLSP, the runs with LSODE-A and LSODE-B were 
required to satisfy the accuracy criteria described in reference 
2. Several of the runs with the other codes also satisfied these 
criteria. For LSODE-B and EPS = 10 ” 4 , two ATOLSP 
values (10* and 10 ^ t: ) satisfied both the accuracy and 
execution time criteria for P2. Table VI includes results 
obtained with both values, to illustrate the effect of ATOLSP 
on the accuracy. However, the error plot given in figure 13 
was generated with ATOLSP — 10 
The maximum percent errors in each species type and 
temperature are given in tables V and VII for PI and P-- For 
each species type the species incurring the maximum error, 
the reaction times at which the maximum errors occurred, and 
the standard solution values for species mole fractions and 
temperature at these times are listed. 

For test problem 1 the runs with EPISODE-A and 
EPISODE-B and EPS>5xlO' 6 predicted little or no 
change in the composition and temperature after an elapsed 
time of 1 ms. Hence, maximum errors of - 100 percent were 
obtained for both intermediates and products (table V). 
Because the temperature and the more active reactants Ht, 
CO, and 0 2 display monotonic behavior (fig. 2) the 
maximum errors in reactants and temperature occurred at 
/. nd , the final time (>l ms) at which the solution was 
generated. Both EPISODE-A and EPISODE-B required only 
six steps to complete the problem, and, because a new step 
size is considered after every successful step, r en j was 
significantly greater than 1 ms. For intermediates and products 
maximum rms errors of 100 percent were obtained at several 
time values; hence, reaction times and standard solution values 
are not given in tables IV and V for these two species types. 
Also, no species name is listed in table V tor either type 
because all intermediate and product species incurred the 
shown maximum errors. 

The solution returned by EPISODE was also found to depend 
on the output stations specified by the user. For example, tor 
some combinations of output times, EPISODE-B (with 
EPS = 1 x 10 6 ) predicted no change in the composition and 
temperature after an elapsed lime of 1 ms. However, by 
stipulating only one output station (at 1 ms), the correct 
solution was obtained. For this problem the run with GCKP84 
and EPS = lxl0~ 2 exhibited serious instability and was 
therefore terminated. For reasons just discussed, no error plots 
for EPISODE-A. EPISODE-B, and the run with GCKP84 
and EPS = 10 2 are presented. 

Similar remarks apply to the results obtained for P2 with 
EPISODE-A and EPS > 5xl0~ 4 , and with EPISODE-B 
and EPS>5xl0 - '\ The run with EPISODE-A and 
EPS = 5xl0~ 4 required only seven steps to complete the 
problem and f cnd was therefore significantly greater than 1 
ms. For this EPS, for exactly the same reasons given for PI , 
the following quantities are not shown in tables VI and VII: 
reaction times at which the maximum rms errors occurred in 
the intermediates and products, intermediate and product 
species incurring the maximum errors, the reaction times at 


which the maximum errors occurred, and the standard solution 
values. 

For EPS = 5x 10 " 4 and 10 3 EPISODE-B successfully 
completed P2 in that correct solutions were returned at t = 1 
ms. However, during heat release they were significantly 
inaccurate. For example, the run with EPS = 5 X 10 4 pre- 
dicted little change from the initial composition and temperature 
until / — 40 /xs when heat release began. In contrast, the 
standard solution shows that heat release is almost over by 
this time (fig. 3). As a result, maximum errors of - 100 percent 
were observed for the products (table VII). This error was 
incurred by several product species at several time steps, 
hence, table VII does not list the product species name, 
reaction time, and standard solution value. Because of the 
difficulties experienced by EPISODE-A and EPISODE-B, 
error plots for EPS > 5xl0~ 4 are not presented. 

As discussed previously, all results with GCKP84 were 
obtained with h 0 = 10 6 s, although its current default value 
is 5 x 10 s. The effects of this change in /t 0 on the accuracy 
and execution time were studied by generating results with 
/ J(> = 5 x 10 8 s. The maximum errors incurred by the 
solutions produced with both /»<> values are given in tables 
VIII to XI. For this study new standard solutions using 
/ )() = 5x 10 s s were established to bias the results in favor 
of the current default value for /t„. Despite this bias, tables 
VIII to XI show that in almost all cases /?„ = 10 6 s produced 
more accurate solutions than h a = 5x10 s. (No results are 
given for PI and EPS = 10 2 because the runs with both A 0 
were terminated due to instability.) For PI the results with 
/i () = 10 f> s were significantly more accurate for all values of 
EPS (tables VIII and IX). Surprisingly, for A„ = 5x 10 * s 
the solution with EPS = 10 4 incurred substantially greater 
errors than those generated with the larger EPS. For P2 the 
differences in errors obtained with the two /j 0 were small tor 
EPS = 10 2 and 10 \ but for EPS = 10 4 , A„ = >0' 6 s 
produced significantly more accurate results than 
/j 0 = 5xl0“ s s (tables X and XI). Finally, the execution 
times required with the two A 0 are comparable for both 
problems and all EPS (tables VIII and X). 

Examination of figures 4 to 17 shows sudden increases in 
the error plots for intermediate species and products. This 
behavior is caused by species reaching values of 0. 1 ppm or 
greater (figs. 2 and 3) and introducing their contributions to 
the rms errors. For example, for PI the intermediate species 
producing the sudden increases in the error plots are H (at 
/ = 2 /is), O (at / = 4 *xs), OH (at i - 4 fts), and N (at t = 
20 /rs). For products the pertinent species are H : 0 (at t - 
2 /is), CO, (at / = 4 /xs), and NO (at / = 15 jis). 

EPISODE-A, EPISODE-B, LSODE-A, LSODE-B, 
GCKP84, and CREK1D all experienced difficulty tracking the 
standard solutions during induction and early heat release when 
the species and temperature change rapidly (figs. 2 and 3). 
The essentially isothermal induction period ends and heat 
release begins when the temperature starts to rapidly increase 
from its initial value. During induction, the reactants and 
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temperature remain fairly constant. EPISODE, LSODE, and 
CREK1D have virtually no errors in the reactants and 
temperature until heat release begins at ~9 jxs (PI) and 3 /as 
(P 2), when the errors in these quantities start to increase. For 
GCKP84, however, the error increases start at earlier times. The 
difference, more noticeable tor PI, is due to the smaller 
reaction times obtained with this code for the onset ot heat 
release. Although, for consistency, we used EPS =10 5 to 
generate standard solutions with GCKP84, values ot EPS as 
small as 10^ 9 (for PI) and 10 ~ 8 (for P2) were found to be 
necessary to achieve tolerance independence of the 
temperature-time trace at early times. 

During induction, the intermediate species and the product 
FEO increase rapidly from negligible initial concentrations. 
The errors in the intermediates and products in this regime 
are, therefore, relatively large. During early heat release 
(r < 15 fis for PI and t < 6 jus for P2), these errors continue 
to remain large as more products are formed and the 
intermediate species continue to change quickly. In this regime, 
the reactants show a sharp decrease, and the temperature rises 
significantly. Many of the ODE’s are unstable (ref. 2), and 
so errors introduced at any step will grow as the integration 
proceeds (refs. 1 and 24). For PI the reactants and temperature 
vary rapidly between 9 and — 15 For P2 the temperature 
rise is not as steep, but the reactants change sharply between 
3 and ~6 /as. Between these times, the errors in the reactants 
and temperature are relatively large (figs. 4 to 7, 9 to 14, and 
17). For P2 the errors incurred at early times are less for 
CREK1D than for LSODE because of the much smaller step 
lengths used by the former code in these regimes (refs. 2, 3, 
and 12). During late heat release and equilibration, however, 
EPISODE, LSODE, and GCKP84 incur much smaller errors. 
In these regimes the ODEN are stable (ref. 2), and so the errors 
decay as the integration proceeds, provided, of course, that 
the numerical method is stable. 

The error plots for EPISODE, LSODE, and GCKP84 
illustrate the dangers of assessing the accuracy of a technique 
(or of a run with a certain value for EPS) by comparing 
solutions at the final time (=1 ms for both problems). Note 
that, although all these codes have negligible errors at the final 
times, the errors can be significant at early times. For example, 
with GCKP84 and EPS = 10 “ 3 the maximum rms error in 
products is over 500 percent for PI . These plots also indicate 
that if the main objective of the calculations is to study postheat 
release phenomena (e.g., NO formation), the use of large error 
tolerances does not result in significant errors. The large errors 
incurred at early times, however, have important implications, 
especially in developing and validating reaction mechanisms. 
A procedure commonly used for this purpose is to compare 
ignition delay times (e.g., time required for the temperature 
to increase by a specified amount) predicted by the mechanism 
with those measured in a shock tube (e.g., refs, 25 and 26). 
The temperature error plots show that caution must be 
exercised in using some of the codes to develop reaction 
mechanisms by applying the above procedure. If, for example. 


we assume that the ignition delay time is the time required for 
a 25 K rise in the temperature, values of - 1 1 and 3.5 fx s are 
obtained for PI and P2, respectively. At these times, the error 
in temperature ranges from 10 to 25 K for EPISODE, 2 to 
5 K for LSODE, 15 to 200 K for GCKP84, and 0 to 10 K 
for CREK1D. 

In contrast to EPISODE, LSODE, GCKP84, and CREK 1 D, 
CHEMEQ incurs virtually no errors during induction and early 
heat release (figs. 7, 8, 15, and 16). Therefore, this code can 
be used to generate accurate ignition delay times. CFIEMEQ 
is superior in these regimes because of the very small step 
lengths that it selects (refs. 2 to 5). However, as pointed out 
by Young and Boris (ref. 10), the continued use of the hybrid 
method used in CHEMEQ results in the global errors 
increasing with time. For example, with CHEMEQ-A and 
EPS = 10 ' 2 , the rms error in reactants has risen to almost 
50 percent for PI (fig. 7) and 25 percent for P2 (fig. 15). The 
situation is worse with CHEMEQ-B (tigs. 8 and 16). During 
equilibration, for CREK ID, also, the errors grow (figs. 9 and 
17) because the formulation used by it in this regime is based 
on that used in CHEMEQ. However, CREK ID incurs smaller 
errors than CHEMEQ for most of the species types and for 
the temperature. For EPS < 10 3 both CHEMEQ and 
CREK ID either are more accurate than or compare favorably 
with LSODE during late heat release and equilibration. 

Figures 4 to 17 and tables IV to VII show the large variations 
in the maximum errors for the different techniques. EPISODE 
and GCKP84 experience the greatest difficulty tracking the 
solutions at early times— rms errors in excess of 100 percent 
are obtained with the two codes. In contrast, the errors incurred 
by LSODE, CHEMEQ, and CREK ID are significantly less. 
Comparisons of the runs with the largest EPS value show that 
LSODE is the most accurate code for PI, and CREK ID tor 
P2. Comparing the errors in the different regimes shows that 
CHEMEQ is the most accurate code during induction and early 
heat release. During late heat release and equilibration, 
however, the other codes are more accurate. 

Examination of figures 4 to 17 and tables IV to VII shows 
that all techniques are tolerance effective in the sense that a 
decrease in the local tolerance generally results in decreased 
global errors. We note, however, that with LSODE not all 
plots show an error decrease with EPS. On the contrary, for 
some runs the error increases with a reduction in EPS (figs. 
4, 5, 12, and 13). This behavior can be explained by examining 
the nature of the error control performed in LSODE. As 
discussed in the section "Computational Procedure, 1 ' the error 
control selected to be performed by LSODE is mixed relative 
and absolute for species mole numbers and pure relative tor 
the temperature. For pure relative error control, the estimated 
local truncation error, d h in species / approximately satisfies 
the inequality 

d, < EPS cr,j (28) 

For pure absolute error control, d, approximately satisfies 



d t < ATOLSP 


(29) 

These two inequalities are approximate because the code 
controls only the rms norm of the estimated local truncation 
errors in all variables and not the estimated local truncation 
error in each variable. 

Equations (28) and (29) show that, since a, « 1 , relative 
error control is more accurate for a given value of the local 
error tolerance. Hence, relative error control is appropriate 
for the two test problems. However, when o i — 0, relative 
error control cannot be used. This problem is resolved by using 
a mixed relative and absolute error control, and 

d, < EPS jot + ATOLSP (30) 

Equation (30) shows that for the error control to be relative, 
ATOLSP must satisfy the inequality 

ATOLSP « EPS a,! (31) 

If ATOLSP » EPS |a, the error control is absolute. In this 
study, we have considered only species with mole fractions 
(x,-) > 0. 1 ppm. This value of X/ corresponds to a, = 3x 10 ~ 9 
and 4x 10~ 9 , respectively, for PI and P2. Hence, for the error 
control to be always relative, ATOLSP must be less than 
3xl0 -9 EPS and 4xl0 -9 EPS, respectively, for PI and P2. 
Only the runs with EPS = 10 ~ 2 for PI satisfy these require- 
ments. Hence, they are the most accurate at early times when 
the mole numbers of many intermediate and product species 
have very small values. Note that for P2 even the standard 
solutions do not satisfy the requirement on ATOLSP. To ensure 
their accuracy, the standard solutions generated with LSODE- 
A and B were checked, respectively, against the solutions 
obtained with LSODE-A and B using EPS = 10 ” 5 and 
ATOLSP = 10 -L \ which satisfy equation (31). These compar- 
isons showed agreement to three significant figures for all 
species with mole fractions > 0. 1 ppm. For LSODE-A agree- 
ment to three significant figures was obtained for all species, 
even those with mole fractions significantly smaller than 
0. 1 ppm. But for LSODE-B the agreement for mole fractions 
< 0. 1 ppm was not good for NO. For all other species, 
however, good agreement was obtained for mole fractions 
> 10~ n . This observation indicates that LSODE-A is more 
accurate than LSODE-B. 

Equation (30) also shows that for given values of EPS and 
ATOLSP, Oj must satisfy the following inequality 

ATOLSP 

^ EP5 = (7min (^ 2 .) 

to achieve relative error control. As a i increases from zero, 
the error control becomes less absolute and is equally relative 
and absolute at a min . For a, > a mjn , the error control becomes 
increasingly relative as a, increases. Hence, the quantity a min 
may be regarded as the value at which, for increasing a h the 


error control starts to change character from being more 
absolute to becoming more relative. Because x, = a,/o,„, the 
value of x, (=x min ) corresponding to a min is given by 

Xmin = ATOLSP/(EPS o m ) (33) 

For PI, x mm = 3xl(T 7 and 3xl0“ 6 for EPS = 10 3 and 
10 -4 and the ATOLSP given in table IV. These values are 
attained by most of the species at t — 5 and 7 /as, respectively 
(fig. 2). Hence, until these times the solutions with EPS = 10 3 
and 10 4 are expected to be worse than or, at best, as accurate 
as the run with EPS = 10 2 . Examination of figures 4 and 
5 shows that the errors in intermediates and products for 
EPS = 10 3 and 10 ~ 4 are worse than those for EPS = 10 2 
until t = 6 and 5 /-is, respectively, for LSODE-A, and until 
t — 7 and 6 /as, respectively, for LSODE-B. In addition, all 
maximum rms and maximum errors, almost all of which occur 
at / > 7/as, exhibit reductions with decreasing EPS (tables IV 
and V). 

For P2 and LSODE-A, x min has values of 2.5x10 \ 
2.5X10 -4 , and 2.5xl0’\ respectively, for EPS = 10 _2 , 
10 \ and 10 -4 . Some of the species never reach these values 
(fig. 3). Hence, the errors do not show much sensitivity to 
changes in EPS (fig. 12 and tables VI and VII). For LSODE- 
B, however, the values used for ATOLSP ensure comparable 
levels of relative error control for EPS = 10~ 2 and 10 -3 ; for 
EPS = 10 “ 4 , the control is more relative in the sense that it 
has a smaller value of ATOLSP/EPS. The errors, therefore, 
display decreases with reductions in EPS (fig. 13 and tables 
VI and VII). The sudden increases in the product errors around 
t = 10 /as were caused by the species NO, which LSODE-B 
had difficulty tracking (table VII). 

The above discussion should be regarded as strictly approx- 
imate because it applies only to the estimated local truncation 
errors, whereas figures 4 and 5 give the estimated global errors, 
which represent the cumulative effects of the local errors. The 
number of integration steps required up to the relevant reaction 
time should therefore also be considered. However, the global 
errors accumulate in a complicated manner from the local errors. 
Other factors that must be taken into account are that LSODE 
controls only the norm of the estimated local errors and that 
different species reach mole fraction values of 0.1 ppm at 
different times. Finally, although we have ignored species with 
x, < 0. 1 ppm, they do incur errors whose magnitudes are 
controlled by ATOLSP and which grow with reaction time in 
the initial combustion regimes, for reasons previously given. It 
is therefore difficult to draw definitive conclusions about the 
ATOLSP to EMAX ratios required for combustion kinetics 
problems. For example, for PI, EPS = 10 " 3 is expected to 
produce more accurate results than EPS = 10 -4 for the inter- 
mediates and products at early times, especially between 5 and 
7 /as (see eq. (33) and the discussion following it), but figures 
4 and 5 show the opposite behavior for both LSODE-A and 
LSODE-B. One conclusion that can, however, be made is that 
care must be exercised in specifying ATOLSP. 


34 


The effect of ATOLSP on solution accuracy is further 
illustrated for P2 and LSODE-B by the results presented in tables 
VI and VII (EPS = 10 ~ 4 ) and XII and XIII (EPS - 10 ). 

Note the significant error reductions obtained by decreasing 
ATOLSP Comparing the errors given for LSODE-B in tables 
VI and VII with those in tables XII and XIII. respectively, shows 
that for the same value (=10“* or 10“ ) of ATOLSP, 

EPS = 10 5 does not produce significantly more accurate 

solutions than the larger EPS. 

Tables XII and XIII show that, although the use of large 
values of ATOLSP can result in significant errors for the 
intermediate species and products, the effect on the temp- 
erature is small. Therefore, if the user is interested only in 
temperature versus time traces at early times, as for examp e 
in developing reaction mechanisms from ignition delay times, 
fairly large ATOLSP values can be assigned. 

The results obtained above indicate that the ATOLSP nee e 
to achieve acceptable accuracy depends as much on the nature 
of the solution as on the value specified for EPS and the mini- 
mum mole fraction to be considered in the error analysis. The 
estimate for ATOLSP given by equation (31) may not be small 
enough, as for example. PI and EPS = 10“ (table IV). On 
the other hand, the estimate may be needlessly conservative^ 
For example, although the intermediate species increase much 
more rapidly for P2 than for PI. larger ATOLSP produce 
results that satisfied the accuracy criteria. Because the value 
needed for ATOLSP is a function of the problem, it can be 
obtained only after the problem is solved. The major problem 
associated with using LSODE to solve chemical kinetic rate 
equations is therefore the trial-and-error procedure necessary 
to obtain the optimal value of ATOLSP, that is, the value that 
minimizes the CPU time while satisfying prescribed accuracy 
requirements. Note that for P2, although the runs with LSODE-B 
and ATOLSP = 10 * and 10' 12 (EPS = 10 “ 4 ) required the 
same CPU time, the latter is significantly more accurate 
(table VI). In contrast, the runs with ATOLSP = 10 , 10 , 

and 10'" required about 2.7, 1.7, and 1.7 s of CPU time 
respectively. The trial-and-error search tor the optimal ATOLSP 
can be time consuming, especially for large systems ot ODE's. 
The use of an extremely small ATOLSP to ensure solution 
reliability can result in excessive CPU times. For examp e 
for P2 the run using LSODE-B with EPS — 10 an 
ATOLSP = l(r M required about 3.4 s of CPU time, in 
contrast, the run with ATOLSP = 10^ 15 required almost 20 s, 
although the solution was not significantly different. 

The error control used in EPISODE and GCKP84 is pure 
relative for species with initially nonzero mole numbers and 
for the temperature; it is. however, pure absolute for species 
with initially zero mole numbers. Since most of the species 
have zero initial mole numbers for both test problems, the error 
control is mostly absolute. Hence, for the same value of EPS, 
EPISODE and GCKP84 are not as accurate as LSODb tor 
solving chemical kinetic rate equations. To achieve comparable 
accuracy, especially at early times when a, is very small 
small values have to be used for EPS (ref. 2). The runs with 


EPISODE and GCKP84 were therefore more expensive than 
the ones with LSODE (refs. 2 and 3). Modifying EPISODE 
to employ the same error control as LSODE produced 
significant reductions in execution times. Preliminary results 
with the revised EPISODE indicate that it is as fast as LSODb. 
For example, for PI the runs with the modified EPISODE-A 
(EPS = 10 2 . 10 "\ and 10 4 and the ATOLSP given in 
table IV) required, respectively, 0.35, 0.41, and 0.61 s. The 
execution times compare very favorably with those requn-ed 
by LSODE- A: 0.37, 0.46, and 0.63 s (refs. 3 and 4). The 
above observations indicate that the error control used in 
EPISODE and GCKP84 is inappropriate for combustion 
kinetics problems. 

Examination of tables IV to VII shows that temperature 
calculation method A does not necessarily produce less accurate 
solutions than method B. On the contrary, for most of the runs 
method A is more accurate for all codes; this result is most 

apparent for CHEMEQ and PL 

To provide a more comprehensive measure for comparing 
the accuracy of the methods examined, we adopted the 
following procedure: For each run, a mean integrated mis 
error, S rms , was defined as 
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1 r'™' 


tend 


M) dt 


(34) 


where f cnd (>lms) is the end of the integration time interval 

and is given by equation (27). 

Equation (34) provides a single quantity that is a measure 
of the average error incurred in solving the complete pro cm. 
The integral in this equation was evaluated numerically using 
Simpson’s rule (e.g., ref. 23), modified for unequal step sizes. 
For runs requiring an odd number ot integration steps, e 
trapezoidal rule (ref. 23) was used on the last two mesh points. 

The effects of h, on S rms for GCKP84 are given in tables 
VIII and X for PI and P2. Except for the run with 
EPS = 10' 3 for P2, ho = 10 ' h s incurred either comparable 
or significantly smaller £ rills than h {) = 5x10 s. Note, 
further, for PI and h Q = 5x 10'* s the substantial increase 
in £ rms when EPS is decreased to 10 (table VIII). 

The variation of 8 rms with ATOLSP is given in table XII 
for P2 using LSODE-B and EPS = 10 " • ™s tab le ,Uus fttfes 
the increasing accuracy obtained by reducing ATOLSP. It al. 
shows that ATOLSP must be chosen carefully, as discussed 

Pr The variations of 8 nns with the user-specified local tolerance. 
EPS are shown in figures 18 and 19 tor an 
respectively. We have included the run with EP 1 SODE-B and 
FPS = 5x 10 4 in figure 19 because it was successfully com 
pL. The 6,„, given in f.gurc .» tor LSODE-B i»d 
EPS = 10 4 was that obtained with ATOLSP 
These figures show that all methods are tolerance effective 
(i c decreasing EPS results in reduced S rms ). For both test 
problems temperature calculation method A is as accurate as 


Mean integrated rms error, % 




method B. In many cases it is significantly more accurate 
especially with CHEMEQ and EPISODE. For P2 and 
EPS =10 , LSODE-B is more accurate than LSODE-A 

because it used a smaller ATOLSP (table VI). 

For the same value of EPS, EPISODE and GCKP84 are 
significantly less accurate than LSODE (figs. 18 and 19) 
because the error control used in the two codes is inappropriate 
for chemical kinetics rate equations. For all techniques, note 
the significant discrepancies between the values specified for 
the user-supplied local tolerance and the errors actually 
incurred. With CHEMEQ-B, a value of EPS = 10 “ 2 (1 per- 
cent) has resulted in an average error of almost 50 percent. 
The relatively large S m]s incurred by CREKID and CHEMEQ 
is due to the difficulties which these codes experienced tracking 
the standard solutions during late heat release and equilibration. 
With LSODE, especially the runs with EPS = 10 -2 the 
correspondence between EPS and S rms is better. These plots 
show that for a given value of EPS, LSODE is the most 
accurate code currently available for solving chemical kinetic 
rate equations. However, for P2, especially with the smallest 
EPS examined, GCKP84, CHEMEQ-A, and CREKID 
compare favorably with LSODE (fig. 19). 

Figures 20 and 21 present the variations of the computational 
work (expressed as the CPU time in seconds) with the mean 
integrated rms error for problems 1 and 2, respectively. Note 
the large differences in the CPU time required by the different 
codes to achieve comparable accuracy. For PI and a Vi percent 
mean integrated global error, the CPU time varies from about 
0.4 s for LSODE-A to over 40 s for CHEMEQ-A. In general 
to produce an order of magnitude reduction in S rms approx- 
imately doubles the computational cost. For both test problems 
LSODE is the most efficient code in the sense that it requires 
the least CPU time to attain a specified accuracy level. 

Figures 20 and 21 show that the CPU times required by 
temperature calculation method A are less than, or compare 



Mean integrated rms error, ^ rms 

Figure -0-— Variation of the CPU time with the mean integrated rms error 
for test problem I . 
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Figure 21. -Variation ol'lhc CPU time with the mean integrated rms error 
tor test problem 2. 


favorably with, those required by method B, This difference 
is most pronounced for CHEMEQ and EPISODE. For 
example, for P2 and a 1 percent 8 nm . CHEMEQ-A required 
only about halfas much CPU time as CHEMEQ-B (fig. 21). 
Note that for EPISODE B the computational work increases 
with increasing error. 

EPISODE-A compares very favorably with LSODE for P2 
(fig. 21). However, the solution generated by EPISODE can 
be strongly dependent on the value selected by the user for 
h 0 and a poor guess can result in incorrect and unstable 
solutions (refs. 2 to 5). It can also result in excessive CPU 
times. For example, the run using EPISODE-A with 
EPS = 10~ 4 and /i 0 = 10 “ s s required about 129 s for P2; in 
contrast, the run with h 0 = 10 7 s required only 0.59 s. 


The accuracy of several codes (EPISODE, LSODE, 
GCKP84, CHEMEQ, and CREK1D) in solving combustion 
kinetic rate equations has been examined in detail. The 
accuracy studies were made by applying the codes to two 
practical’ combustion kinetics problems. Both problems described 
adiabatic, homogeneous, gas-phase chemical reactions and 
included all three combustion regimes: induction, heat release, 
and equilibration. 

During induction and early heat release, when the species 
mole numbers and temperature change rapidly, EPISODE, 
LSODE, GCKP84, and CREK1D had difficulty tracking the 
solutions. The errors incurred by EPISODE and GCKP84 in 
these regimes were significantly larger than those incurred by 
LSODE and CREK1D. In contrast, the solutions generated 
with CHEMEQ displayed virtually no errors during induction 
and early heat release. However, during late heat release and 
equilibration, the errors obtained with CHEMEQ increased 
significantly. In these regimes, the other codes were more 
accurate. 

Among the codes examined, LSODE was the most accurate 
for solving chemical kinetics problems. This study has also 
shown that LSODE is the most efficient code, that is. it 
required the least execution time to attain a specified accuracy. 
The major difficulty associated with its use is the trial-and- 
error procedure necessary to obtain optimal values for the local 
absolute error tolerances tor the variables. A poor guess for 
the absolute error tolerance can result in excessive execution 
times or in seriously inaccurate solutions. 

An important conclusion is that calculation of the temperature 
by solving the algebraic enthalpy conservation equation can 
be more accurate and efficient than integrating its differential 
equation. 
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Appendix A 

Description of Codes Studied 


The ordinary differential equations (ODE’s) (2), (9) and (12) 
describing homogeneous gas phase chemical reactions can be 
generalized as follows: 


y, = ^T=f((y k )) i.k= l n 

dt 

ydb) = given 


(Al) 


where for temperature calculation method A (see the section 
“Evaluation of Temperature’’) 


= a, i = 1 N s 

N = N S 

for temperature calculation method B 

= a, i = 1 ,...,N S 
yN s + 1 = t 
N = N s + 1 

and for the code GCKP84 


(A2) 


(A3) 


y t = Oi i = 

> J N S + 1 := V >W A .+ 2 = P yN s + 3 = T 
N = N s + 3 


(A4) 


In vector notation equation (Al) becomes 


y = 



y(*o) = given 


(A5) 


where the underscore is used to denote a vector quantity. A 
matrix is denoted by a boldface letter. This notation is used 
throughout this appendix. In equation (A5) the /V-dimensional 
column vectors y and £ contain the dependent variables and 
their temporal derivatives, respectively. 

The initial-value problem is to determine values for (y,j at 
one or more times in a prescribed integration interval, given 
[fi] and the values [y^ ( t 0 )} at the initial time t 0 . We now 
describe the codes studied in the present work and how they 
solve the above problem. 


EPISODE AND LSODE 

Both these codes use linear multistep methods of the form 
(refs. 6 to 9) 


*1 *2 
i j= o 

(A6) 

where Y in is an approximation to the exact solution >>,(*„), 
f, n (equal to fi([Y kn })) is an approximation to the exact 
derivative y { (t n ) (equal to/((y* (/„)})), and the {a J>n \ and 
(&o, n > 0) are associated with the particular formula selected 
by the user. The options include a variable-step, variable-order 
implicit Adams method (suitable for nonstiff problems) of 
orders 1 to 12, and a variable-step, variable-order backward 
differentiation formula (BDF) method (suitable for stiff 
problems) of orders 1 to 5. As discussed in the section 
“Methods and Codes Examined,” the BDF method was more 
efficient for the problems examined in this study. Therefore, 
the discussion is restricted to this method. For a BDF method 
of order q , K x = q, K 2 = 0, and equation (A6) reduces to 


<7 

^i.n ~ ^ + ^n0O,nfi,n 1 = \ • ,N (A7) 

j= I 

The step length h n can vary from one step to the next in 
EPISODE but is held constant for q + 1 consecutive successful 
steps in LSODE. Hence, for EPISODE, {ctjJ and (/3,-J can 
vary from one step to the next, but in LSODE they are 
predetermined constants corresponding to the order used. 

Both codes use a predictor-corrector process to solve for 
Y„. An explicit method generates a predicted value, Fj 01 , 
which is then corrected by iterating equation (A7) to convergence, 
that is, the improved estimates Yl m] (m = are produced 

until the magnitude of the difference (Fj" ,] - ,J ) T in 
EPISODE, or (h n P n m ^ - h n Y^~ 1] ), in LSODE, approaches 
zero within a specified accuracy. Here, Y}^ and Y are, 
respectively, the approximations generated for Y n and^ on 
the m th iteration, the integer M is the number of iterations 
required for convergence, and Yj, m] is accepted as the numerical 
solution at t n , provided it satisfies a prescribed local accuracy 
requirement. At each iteration m, h n Y}™ ] is computed in LSODE 
from Y}^ via the relation 

Ylr ] = Y l <Xj.nY n ,j + h n (3 0 ' n Ylr ] (A8) 

j= 1 

so that the pair (T^, h n Y)™^) satisfies the BDF method (eq. 
(A7)) exactly. The predicted values of Y n and denoted 
by h n Y^\ also satisfy equation (A8). 

The predicted quantities yj 0] and h n Y® ] are obtained by a 
<? th -order Taylor series expansion as follows: The history of 
the solution is maintained in the Nordsieck array (which is 
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a Taylor series array) z„ of size N x (q + 1) (e.g., ref. 1). 
The / th row (7 = l,..., AO contains the q + 1 elements Y in , 

h n Y lfr h;J2\ Y lf JYJJq! Yj%\ where Y$ is the approximation 

to (cl 1 Y, /clt J ), . The (q + 1) columns of z„ are numbered from 

0 to q , and the 7 th column (j = QA q) y w-hich will be 

denoted by the vector z„(j)< contains the vector h„Y^ ] lj\ of the 
/ h -ordcr scaled derivatives. If z, r _, has been obtained, the 
predicted history matrix, zj, ()1 , at t n is given by (ref. 1) 

z, 1 , 01 = z„_,A(< 7 ) (A9) 


where A(q) is a (q+ 1) x (</+l) matrix, with element A fk {q) 
given by 


Aj k (<i) = 


0 7 < O 

(i) j*k) 


jrk = 0,1 q 


The binomial coefficient, (j.), is defined as 


and LSODE use different calculation procedures, the [I iw „\ 
values are, in general, different in the two codes. For 
EPISODE, I„(q) depends on the method order and the step 
length history, satisfies / 0 ,,(</) = 1 and /| „ = l//3 0 ,„ and has 
to be recomputed at the start of each step. For LSODE, l n (q) 
is a function of only q y satisfies l^„(q) = (3 0n and /,„(</) = 1 , 
and has to be recomputed only when the method order is 
changed. 

To correct the initial estimate J^ 01 (i.e., to solve equations 
(A7)), both codes include a variety of iteration techniques. 
For combustion kinetics problems the most efficient is the 
Newton-Raphson iteration (ref. 2), which is given by the 
recursive relation 


p(yd;"+" _ 


"<\ 


£ + h,A J ( H" 1 ) - Vi 

./=! 

(A 12) 


for m = 0,1 A/— 1. The NxN iteration matrix P is 

given by 


(7 ) = lL 

k! (j — k)! 

Thus, the predicted array z, 1 , 01 is obtained by a simple q' h - 
order Taylor series expansion by using equation (A9). The 
matrix zj, 0 * contains predicted values of Y n and its scaled 
derivatives up to order q , the current method order. Note, 
however, that because a c/ lh -order Taylor series expansion 
method is used, zj^iq) = z n -\(q). 

The estimates Y)^ and, in LSODE, h n Y}^ (m = 1 
are generated, as described below, until the iteration 
converges. The local error test is then applied and, if passed, 
the Nordsieck history matrix z„ is constructed by using the 
relation 

■/,„ ~ zj, 0 ' + e„l „(q) (A 10) 

where 

_ y\M\ _ ydO) 

L ti 1 n } n 

in EPISODE, and 

e - h - li T 101 

*-n Il n A n u n 1 n 

in LSODE. The (q + 1 ) th -dimensional vector l n (q) 

I'M) - (/()., ,(</). (All) 

contains the method coefficients for the Nordsieck history 
formulation of the c/ h -order BDF method. Because EPISODE 


P = I - h„0 {)jt S (A 1 3) 

where I is the identity matrix and J is the Jacobian matrix, 
w'ith element given by 

J'j = df/dy, i’j = 1 A 

For this method, much computation time is required to form 
the Jacobian matrix and to perform the linear algebra necessary 
to solve equation (A 12). To reduce this computational work, 
P is not updated at every iteration. For further savings, it is 
updated only when it has been determined to be absolutely 
necessary for convergence. Hence, the iteration matrix is only 
accurate enough for the iteration to converge, and the codes 
may use the same matrix over several steps of the integration. 
In any case, both EPISODE and LSODE update P at least 
every 20th step. The linear algebra required to solve equation 
(A12) is performed by the LU method (e.g., ref. 27), rather 
than by explicitly inverting the matrix, which requires 
prohibitive amounts of computer time (ref. 23). 

Convergence of the estimates is ascertained as discussed 
below. EPISODE constructs a vector F max , which depends 
on the user-specified value for the local error control IERROR 
as follows: 

IERROR = 1 (absolute error control): 

^,ax,= l ' = 1 N 

IERROR = 2 (pure relative error control): 

Y nm , ■= IX,,-. I ' = 1 N 



IERROR = 3 (semirelative error control): 

I'max, = max! \Y in _ 1 1, I^.„_ 2 |l for i that satisfy y,(/ 0 ) * 0 

- maxjl, |y /)/( _||j for / that satisfy >v(r 0 ) = 0 

(A 14) 

The test for iteration convergence is based on the successive 
differences (Yl m] - L*" 1-11 ), as compared with T max and the 
user-supplied local error tolerance parameter EPS. Conver- 
gence is said to occur if 


where C m ( = b„,/b m _ \ ) is the convergence rate. This assum- 
ption is used to modify the convergence tests (eqs. (A 15) and 
(A 16)) as follows: 

for EPISODE 

K < C E EPS 

where 

b m S m min (1, C m ) 


b 


rn 



(A 15) 


C;„ = max (0.1 C„,_|, C,„) 


and for LSODE 


In LSODE an error weight vector EWT is constructed as 
follows: 


EWT, - RTOL, \Y it , | + ATOL, i = 1 A/ 

where RTOL, and ATOL, are, respectively, the user-supplied 
local relative and local absolute error tolerances for the / ,h 
component. Both RTOL and ATOL can be specified either 
as a scalar or an array, as discussed in the section “Computational 
Procedure/’ The value of the user-supplied parameter ITOL 
indicates whether RTOL and ATOL are scalars or arrays. 
ITOL has four possible values which correspond to the types 
of RTOL and ATOL as follows: 

ITOL = 1 : scalar RTOL and scalar ATOL 
ITOL = 2: scalar RTOL and array ATOL 
ITOL = 3: array RTOL and scalar ATOL 
ITOL = 4: array RTOL and array ATOL 
The convergence test is based on the successive differences 
(h n Y}™^ - h n Y} } m ~^) as compared with EWT , and is given by 


b 


m 



/, y.M _ /, y> 

fl n 1 i,n ri n I t,n 

EWT, 



(A16) 


The factors C E and C L in equations (A 15) and (A 16) are 
chosen to make the convergence tests consistent with the local 
truncation error tests. In particular, C E = 0. 1 3 E (q) and 
C L = 3 L (q)!2(q -I- 2), where 3 E (q) and 3 L (q) are the test 
constants used, respectively, in EPISODE and LSODE for the 
local error test (eqs. (A27) and (A28)) and where the variable 
q indicates the method order. 

If convergence is not achieved after the first iteration, the 
codes anticipate the magnitude of b m one iteration in advance 
by assuming that the estimates converge linearly. Thus, 6 m+) , 
which does not yet exist, is estimated by 


= b„ 


X-I 




- Q 

where 

b t ' u = b m min (1, 1.5 C' t ) 

C; n - max (0.2C ffl _ ,, C m ) 

Now, at least two iterations are required to compute C m . For 
the first iteration, C,J, is set equal to 1 in EPISODE and equal 
to the last value of C m from the previous step in LSODE. For 
the first iteration of the first step and after every update of 
the Jacobian matrix, LSODE sets C ^ equal to 0.7. 

If the corrector iteration fails to converge in three iterations, 
h n is reduced by a factor of four if P is current and the step 
is retried; otherwise, P is updated at y = Lj, 01 , and the step 
is retried. The same corrective actions are taken by LSODE 
if C m > 2 after the second iteration. In the event of a singular 
iteration matrix, both codes reduce h n by a factor of four and 
attempt the solution with the new step length. The integration 
is abandoned if either the step size is reduced below a minimum 
value (both codes) or 10 convergence failures have occurred 
(LSODE). 

If the corrector converges after M (<3) iterations, an 
estimate of the local truncation error is made, as described 
below. For a BDF method of any order k , the local truncation 
error, d n (k) at t„ is given by 

d„(k) = e k+l h k ;' y«+"(t n ) (A 17) 

where the variable k denotes the method order and the constant 
C* +1 depends on the method formulation. For the variable- 
step method used in EPISODE (ref. 6), 


y= i 

(*+!)//,.„(*) 


(A 18) 


40 



where 


and for LSODE 


= t„ — k 

h„ 


, , , q'k,Aq) ki,M) 

d„{q) = — - — — e„ = e, 

q+\ q+\ 


and /, „(k) is the second element of the ( k + l) lh -dimensional 
coefficient vector l„(k ) for the £ lh -order method (see 
eq. (All)). For the formulation used in LSODE (ref. 28) 


because q\l qjx {q) = / 0 .„(<7) for the formulation used in 
LSODE (ref. 28). 

The local error tests used in the two codes are as follows: 
For EPISODE 


The error d n (q) in the </ ,h -order method (i.e., k = q) used 
on the / 7 th step (i.e., (/„ _ i • /,, ) ) is estimated as follows: As 
discussed previously, the last column of z„, z n (q) contains the 
vector h l ,[Y^ ] lq\ and that of zj 01 , z^Hq) contains the vector 
| iq\. The difference of z n (q) and gives 


Z n {q) -z}?'(q) Y^ {) 4- 0(hy 2 ) (A21) 

q- 

by using the mean value theorem for derivatives. From 
equation (A 10) the above difference is seen to be equal to 
l qn (q) e n , which, upon substitution into equation (A21), gives 

/>r'^ +l) - q! q)<?„ (A22) 

if higher-order terms are neglected. This approximation is used 
in LSODE. EPISODE, however, takes into account errors in 
the past values and uses the following expression (ref. 6): 


" / d in Y\ l/2 ? 


< EPS 


which, upon using equation (A25) can be written as 


s Eir 

1 j = | \ ^ max 


3 t: (q) EPS 


where the test constant 3 E (q) is given by 


3 E {q) =i\,(q) 1 + n (fz7 


j- 2 V" n -.A 


And for LSODE 


hr'Y^'^-e, 
y „ 


i y f^-Y X 1 9 

N “ \E WT 7 j 


(<?+>)/ 


+ n ; 


-j 


n — 1 1 n - /, 


and where is given by equation (A 19). 

By substituting the above expressions for h tl Y f \ q ^ ]) and the 
appropriate equation (A 18) or (A20) into equation (A 17) (with 
k = q) and simplifying the resulting expressions, we obtain 
the following estimates for d n ( q ) : 

For EPISODE 


By using equation (A26), the above inequality can be expressed as 


i N / \A I/2 

- Y ( \ 

N U \EWT,/ J 


where the test constant 3 L (q) is given by 


/,.»(?) 


n — 

f = \ \<n - 1 - t„-j 


e„ (A25) 


3 ,.(q) = 


htjAq) 
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If the error test fails, the following corrective actions are 
taken. In EPISODE h n is reduced so that equation (A27) is 
satisfied (see eqs. (A29) and (A30)), and the step is retried. 
If the results with the new step length do not pass the error 
test, h n is reduced by a factor of five. After three and more 
error test failures EPISODE reduces the step length by a factor 
of 10 and reduces the method order by one if it is greater than 
one. If an error test failure occurs with q — 1 , the Nordsieck 
history matrix z, ( _i is reconstructed from Y n _ x and /( Y n _ j ) . 
After the first error test failure, LSODE reduces h n and/or 
q by one and then retries the step. If the error test is again 
not satisfied, h n is reduced by a factor of five. After three and 
more such failures, the method order is reduced to one if it 
is greater than one, the step size is reduced by a factor of 10, 
and the step is retried with a new Nordsieck history matrix 
z„_|, which is constructed from Y n _ { and f(Y n _ x ). Both 
codes abandon the integration if h n is reduced below a 
minimum value. The maximum number of error test failures 
allowed is seven in EPISODE and 10 in LSODE, after which 
an error exit is taken. 

If the error test is passed, the step is accepted as successful, 
and the entire Nordsieck history array. z„, at t n is updated by 
using equation (A 10). 

Periodically, both codes attempt to change the step length 
and/or the method order to minimize computational work while 
maintaining prescribed accuracy. After every step on which 
no convergence test or local error test failure occurs, 
EPISODE attempts to use a larger step length at the same 
method order. The new step size h'(q), where the variable 
q denotes the order to be used on the next step, is chosen such 
that it exactly satisfies the local error bound (eq. (A27)) by 
assuming that the highest derivative remains constant. Then, 
because d n varies as h% +l (eq. (A17)), 


maximum order, </ max , the choice q + 1 is not considered. 
Also, if q = 1 , the choice ^ — 1 is rejected. For each method 
order q ' the step size h'(q') that will exactly satisfy the local 
error bound is obtained by using the procedure outlined above 
for q' = q (eq. (A29)). 

For the case q' = q — 1, d n (q— 1) varies as h^y ^ q) (r„) (eq. 
(A 17)), which is equal to q\z n (q). The local error test for 
q' = q — 1 is as follows: 

For EPISODE 



3 £ (<7-l) EPS 


where Zj,„(q) is the / th element of z n (q) and (ref 6) 


3,(?-l) 


i\(q-i) 

q-\ 

FL 

j= i 


And for LSODE 


A,-, 


i N /, ( \ 1/2 

N £ V EWT, / ] ? 

< 1 

3i(?-l) 




h’ (q) 


q+\ 


(A29) 


where r is the ratio of the step length to be attempted on the 
next step to its current value and the subscript “same” 
indicates that the current order ( q ) is to be used on the next 
step. To allow for inaccuracies in the error estimate, certain 
safety factors are built into the calculation procedure for h'(q) 
to produce a smaller value than that given by equation (A29). 
The formula used in EPISODE for r sa me is 


1 

T;ame — — - 

(5 D q ) q+ ' 4- 10~ 6 


(A30) 


To increase the efficiency, both codes consider changing 
the method order to q — 1 or q -V 1 at periodic intervals. After 
an unsuccessful step or when the current order equals the 


where (ref. 28) 


\(q~ 1) 


1 

(</- 1).' 


The step length ratio, if the order is to be reduced to q — 1 , 
is then given by 




h'(q-\) 



(A31) 


where the subscript “down” indicates that the order is to be 
decreased. If q= 1, r down is set equal to zero because q 
cannot be decreased. 

For the case q' = q + 1, d n (q+ 1) varies as h% +2 y iq + 2) (t n ), 
which is estimated by differencing the quantity h q + ] Y i<i+i) 
over the last two steps and then using the mean value theorem 
for derivatives. For EPISODE equation (A23) gives 
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e n ( h tl \ I ^ fo(f+ I y-(</+ I ) _ /^/+ 1 y(<y+ 1 ) 

7 m \^M-i/ 7 m- i 

= hr 2 rt q+2) + o(*r 3 ) 

where y„ is given by equation (A24). For LSODE, equations 
(A22) and (A26) show that the approximation for h? } + ] Y , \ q + 1 } 
is given by „ ( q)e tv Because the methods used in this code 
are based on a constant step size, the quantity l 0n (q)[e n — |] 


where the subscript "up'’ indicates that the order is to be 
increased. If q — q mdX and in LSODE after a failed step, r up is 
set equal to zero to prevent an order increase. 

For reasons given previously certain safety factors are built 
into the step length ratios (eqs. (A29), (A31), and (A32)). The 
formula used in EPISODE for r^ mc is given by equation (A30); 
the other two ratios are computed as follows: 


/o„(<7)k, - = hr 2 y,\ q+2) + o(*r 3 ) 


(5D q _ | ) q + IQ’ 


The local error test for q ' = q + 1 is given by 
For EPISODE 


(10 o v+l r ‘ + io 


M L w 


The formulas used in LSODE to calculate the step length ratios 


where (ref. 6) 


3 £ (</+l) EPS 


1.3 (£>„-,)*+ IO' 


3 £ ( 9 +l )= ii±2)M9±l> 

£</ + I j = 2 ~~ 1 *n-j. 


1.2 j(D (/ ) </+1 + 1(T 6 


And for LSODE 


where (ref. 28) 


\uh\ ewt - / / ? 

D q+l ^ -1 ^ 

3t(9+l) 


3 t (?+l> = 




The step length ratio, if the order is to be increased to q + 1 , 
is then given by 


1.2 (D„ +l ) ^ + 10- 


The order corresponding to the maximum step length ratio 
r = max(r down , r^j^) and the step length ratio r are selected 
to be attempted on the next step if, after a successful step, r > 1.3 
(EPISODE) or LI (LSODE); otherwise, both changes are 
rejected. After a failed step, the order is decreased in LSODE 
if r down > r^ mc ; however, r is set equal to one if it is greater 
than one. Several additional tests are performed on r before the 
step length to be attempted next is selected. These tests may be 
summarized as follows: 

For EPISODE 


h'(q+D ( 1 y + 


\ ^llktx ( h mm 

r - min 1 , r lnax , max r, — — . r min 

( h„ V h„ 
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where the arrow denotes the replacement operator. 
And for LSODE 


r — min 



''max' ™ax 



In EPISODE, /t min and /i max are set equal to, respectively, h (h 
the user-supplied value for the step length to be attempted on 
the first step, and 10 (f out - t oulM ) where r out is the current 
time at which the solution is required and r outold is the 
previous value of r out . On the first call to EPISODE r oul is set 
equal to f 0 , the initial value of t. In LSODE, however, h mm 
(default value = 0 ) and /? max (default value = oo) are user- 
supplied optional input parameters. The quantity r min (used 
only in EPISODE) is set equal to 0. 1 , and r max depends on 
the code. In EPISODE, r max is set equal to 10 for the first 
10 integration steps; thereafter, it is set equal to 1.5. In 
LSODE, r max is normally set equal to 10; for the first step 
length increase following either a convergence or local error 
test failure, it is set equal to two. In both codes for the first 
step length increase for the problem r max is set equal to 10 4 
if no convergence or error test failure has occurred. 

After the step length ratio r has been computed, the step 
length h' to be attempted on the next step is given by 


adds a column of zeros, representing z n (q + 1 ), to z„. 
LSODE sets z n {q + l) equal to h\ ,7 +l T„ U/+I) % +1)!, which 
is equal to l (j Aq)e n /(q + 1 ) (see eq. (A22)). Both codes then 
rescale all q + 2 columns of z„ by powers of h ' I h„ to account 
for any change in the step size. 

The solution values at prescribed output times r olil ,, / OUI 
are obtained quite easily from the history array. For each 
output station f nul , the codes continue the integration until the 
first mesh point n for which t n > r oul and then compute the 
solution at t ou( by a 1 th -order Taylor series expansion 
about t n : 


yu out) 


(t - t V q ’ 1 + ] 

L^hL-_ n) yjj) _ ^ 


j = o 


7=0 



(A33) 


where q t ' J+l and h' +l are, respectively, the method order and 
step size to be attempted on the next step. 

Both codes start the integration with a single-step, first-order 
method because information is available at only the initial 
point, / 0 . The Nordsieck history matrix z 0 at t 0 is constructed 
from the initial conditions y(/ 0 ) and the ODE's as follows: 


Zo(0) — Y 0 — y (/o); so( O = — hofiYo) 


h ' - rh }V 

Changes in method order (and step length in LSODE) are 
attempted only after S successful steps with the same order (and 
step length in LSODE), where S is normally set equal to q + 1 . 
However, if an unsuccessful step occurs, the step length and/or 
order may be reduced. Following a failed error test or a failed 
convergence test, if P is current, EPISODE resets S equal to 
2 if it is less than two, but LSODE resets S equal to q + 1 
irrespective of its current value. If three or more error test failures 
occur on any one step S is set equal to five in LSODE and either 
q + 1 (it q > 1 ) or 10 (if q = 1) in EPISODE. Following a step 
for which the method order is not changed EPISODE sets S equal 
to 2. If method order and step length changes are rejected because 
r < 1.1, LSODE sets S equal to 3. 

After every S - 1 successful steps, if q < <y max , EPISODE 
saves e and 7 , and LSODE saves e„ so that r up can be 
computed. To minimize storage requirements, the vector e is 
saved as the (<? max + l) lh column of z. 

If the step size and/or method order is changed on the /z th 
step, z„ has to be modified. For the case q ' = qji' ^ h„, the 
7 th column (j = 0,1,.., < 7 ) is multiplied by (h'lh n ) J . For the 
case q' = q - 1 , h' ^ h„, the last column of the old z n is 
ignored because it is not needed on subsequent steps. In 
addition, EPISODE adjusts the first q columns to reflect the 
reduced set of data represented by z„ (ref. 6 ). In both codes 
the above scaling by powers of (h'/h n ) is performed on the 
first q columns. For the case q f = q 4 - 1 , h' ^ h n , EPISODE 


where h 0 , the step size to be attempted on the first step, has 
to be supplied by the user to EPISODE. In LSODE, however, 
/? 0 is an optional input variable and is computed by the code, 
unless the user has specified a value for it. 

GCKP84 

GCKP84 is a general-purpose chemical kinetics code 
designed to solve a wide variety of problems (ref. 15). It uses 
the integration technique developed by Zeleznik and McBride 
(ref. 18). As implemented in GCKP84 the integration 
algorithm is an extensively modified version of the GEAR 
package (ref. 19), which is similar to LSODE. In particular, 
GEAR includes the two linear multistep methods discussed 
previously. The methods are based on a constant step length, 
and the method coefficients [l J (eq. (All)) have the same 
values as in LSODE. Hence, l n ( q) is a function of only the 
current method order q , satisfies l 0n (q) = (3 0n (see eq. (A7)) 
arj d h, n (q) = 1 , and has to be recomputed only when the 
method order is changed. GCKP84 uses the same two linear 
multistep methods but the maximum method order is different: 

1 1 for the implicit Adams method and 8 for the BDF method. 
The methods are also implemented differently as discussed 
below. For reasons given previously we restrict discussion to 
the BDF method (eq. (A7)). 

As in EPISODE and LSODE, GCKP84 maintains the 
solution history in the form of the Nordsieck history array, 
z. The array z„ at the current time t n is obtained by using a 
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predictor-corrector process. The prediction step is performed 
in two stages. First, an initial estimate for z, 1 , 01 is computed 
via equation (A9); that is, the result of the prediction step used 
in LSODE and EPISODE (and GEAR) serves only as an initial 
estimate for z'° ! in GCKP84. Second, the above result is 
modified by means of an expression similar to equation (A 10), 
as follows: The difference (z„( 1) - 4 W, <D) (equal to e„ in 
GEAR because l x [q) = 1, eq. (A 10)), that is. 


P is used for a maximum number of 20 steps. At each iteration 
the approximation h„Y}!^ to h ti f n is computed by using equation 
(A8). If any F, 1 "' 1 < 0, the iteration is abandoned. The step 
length is then reduced by a factor of two, and the step is retried. 

Convergence of the estimates is said to occur if any of the 
following three tests, which are applied in the order they are 
given here, is satisfied. The first test involves the magnitude 
of the successive differences ( T’/"’ - T)/" - ' 1 ): 


= hX - A„fi 01 (A34) 

may be regarded as the error in the (/'''-order predictor relative 
to the converged array z„. Equation (A 10) gives the history 
matrix z„ by adding the remainder term associated with using 
the c/ h -order predictor. Of course, since e n can be computed 
only after the converged solution is produced at t fr the above 
procedure cannot be used. However, since c n - \ is available, 
it can be used to improve the initial estimate given by equation 
(A9). However, for additional accuracy improvement, 
GCKP84 uses the quantity E obtained by accumulating the 
errors [e’J (see eq. (A35)). The quantity E n may be regarded 
as the estimated global error in /? /f Fj, 01 . Since E n is not known 
at the start of the step, GCKP84 uses to improve the 
estimate given by equation (A9) as follows: 



where Y tmxJ is given by the expression used in EPISODE for 
semirelative error control (eq. (A 14). The second test is based 
on the size of the current estimate for e n relative to the size 
of the current estimate for E n (see eqs. (A34) and (A35)): 


kyi: a - 

hj}::' - un," + ('f) 

\fh:/ 




- z! n| + 



E n -\l„(q) 


If for any / the denominator in the above summation is less 
than 10 so , it is set equal to 10 2 . The third criterion is 
based on how rapidly the iteration is improving the solution 
and is given by 


where h h , which is normally equal to h n _ j , is the step size 
that E n _ i is based on and the term (/?„//?/. ) ^ 1 accounts for 
this fact: the exponent q + 1 arises because the current order 
is q and the local error varies as 1 (see eq. (A 17)). On 
the first step, E fJ _ x (equal to E 0 ) is set equal to zero because 
Y {) (equal to f(y(t {) ))) is known exactly. 

After the prediction process is performed, the code checks 
the | Y}"} j for negative values. Because it is physically impossible 
for species concentrations, temperature, density, or velocity 
to be less than zero, the results of the predictor step are rejected 
if any Y$ < 0. Also, for each variable i for which the above 
condition is obtained, E iM _\ is reset to zero it it is less than 
zero. The step length is then reduced by a factor of two and 
a new z),° ! is generated. The above procedure is repeated until 
either all predicted solution components are nonnegative or 
the step length is reduced below a minimum value, /? min , in 
which case an error exit is made. 

To correct the initial estimate GCKP84 includes a variety 
of iteration techniques. For reasons given previously the 
discussion is restricted to the Ncwton-Raphson method. The 
procedure used to generate the improved estimates Fj/" 1 
(m = 1,2,...) is exactly the same as that described for LSODE: 
solve equation (A12). The iteration matrix P (eq. (A 1 3)) is 
only accurate enough to achieve convergence, but the same 


| | 

which can be applied only after tw'o iterations. However, the 
third test is applied only after five iterations and that too only 

if S, H — 5. 

If convergence is not achieved after four iterations, the 
iteration matrix P is updated at v = and the correction 
process is retried. This procedure is repeated four times, after 
which, if the estimates have not converged, the step length 
is reduced by a factor of two and the step is retried. The same 
corrective actions are taken if on the Fifth or subsequent 
iteration 6 m > 6,„- X . The above cycle of updating P every 
four iterations and then reducing h n by a factor of two after 
four such updates is repeated until either convergence is 
obtained or the step length is reduced below /i mjn , in which 
case an error exit is taken. 

After corrector convergence the local error test is applied. 
This test is based on E tr which is estimated by using 
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and can be written as 




EPS 3 G (q) 



EPS 3 G (q) 


(A 36) 



EPS 3 c (fl+l) 


where EPS is the user-supplied local error tolerance and 
3 G (q) (equal to 2il 0n (q)) is the local error test coefficient for 
order q. 

If the error test fails, the error vector j is updated by 
using equation (A35), and h E is set equal to h n because E n -\ 
is now based on h tr The code GCKP84 then reduces the step 
size and/or the method order by one and retries the step. After 
three and more error test failures, the method order is reduced 
to one if it is greater than one, and the step length is set equal 
to A new Nordsieck history matrix at t n _ | is constructed 
from K w _ j and E n - \ set equal to zero, and the 

step is retried. After seven such failures or if h n is reduced 
below h m in , the integration is abandoned and an error exit is 
made. 

If the error test passes, the step is accepted as successful, 
the Nordsieck history array z„ is updated by using equation 
(A 10), E„ is computed by means of equation (A35), and h E 
is set equal to h n . 

To increase the efficiency of the integration, the code 
periodically considers changing the method order to q - 1 or 
q + I . Of course, if q = 1 , the choice ^ — 1 is not considered. 
After an unsuccessful step or if either q is equal to the 
maximum method order, c/ max , or D q > 4 t G {q), the choice 
q - h 1 is rejected. For each method order q' the step size 
h'(q') is computed from an estimate of the local error in a 
manner similar to the procedures used in EPISODE and 
LSODE (eqs. (A30) to (A32)). For each method order q' 
GCKP84 computes the step length ratio r(q') as follows: 


r(q’) = 


h'(q') 

h„ 


(A3 7) 


~{D a D.) q ^ 1 + 10~ 6 
3 q 


where 


1 jy (Zi.n(q) ~ 0>5l q .n(q)Ei,n-\ XlX] ~ 


L /V r 


D n 


i=l 


EPS 3 g ( <7 — 1 ) 


and 


The local error test coefficients 3 G (q— 1) and 3 (; (q+ 1) for 
orders q— 1 and q+ 1, respectively, are given by 

3c(9-n =2 


and 


3 C (?+1) — 2 (q + 2)/l 0n (q) 


The quantity D, in equation (A37) is set equal to 10, unless 
A (I _, > 10 25 /vW, in which case it is set equal to 0. 1 . If 
A„ is also greater than 10 ~ 25 /y/N, D- is set equal to A,,/ 
A„_|. Finally, if D : is less than (0.25)' /+1 , it is set equal to 
this quantity. 

The order corresponding to the maximum step length ratio 
r = max(r(q — 1), r(q ), r(q+ 1)) and the step length ratio rare 
selected to be attempted on the next step if /* > 1.1 after a 
successful step; otherwise, both changes are rejected. After 
a failed step, q is decreased if r(q— 1) > r(q); however, r is 
set equal to 1 if it is greater than 1 . The following additional 
tests are performed on r before the step length h' (equal to 
rh n ) to be attempted next is selected: 



The minimum, /t min , and maximum, /* max , step sizes are, 
respectively, set equal to /i 0 , the user-supplied value for the 
step length to be attempted on the first step, and 10(r (Xlt — 
Wold)- 0° ^e first call / out is set equal to f 0 - The quantity 
r max is set equal to 10. For the first step length increase 
following either a failed convergence test or a failed error test, 
it is set equal to two. However, after three or more error test 
failures, it is set equal to min (10 4 , h n /h mi n ), thereby ensuring 
that the new step length equals h mm . For the first step length 
increase for the problem, r max is set equal to 10 4 if no con- 
vergence or error test failure has occurred. 

Changes in method order and step length are attempted only 
after S successful steps with the same order and step length, 
where S is normally set equal to q + 2. However, if an 
unsuccessful step occurs or if D q > 4 3 c (g), the step length 
may be changed, and the method order may be reduced. 
Following a failed convergence or local error test, 5 is set equal 
to q + 2. After three and more error test failures, S is set equal 
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to three. If method order and step size changes arc rejected 
because r < 1 . 1 , S is set equal to 10. Finally, the successful 
step counter is increased by one only if convergence is obtained 
in eight or fewer iterations. 

If the step size and/or the method order are changed on the 
/i ,h step, z„ has to be modified. For the cases q' — q - 1 and 
q' = q, the modifications are made exactly as in LSODE 
(described previously). For the case q' = q 4- 1, z n is first 
augmented by a column containing the vector l qM (q ) 
E n / {q + 1), which is approximately equal to h% + l Y,\ q+ 1 ) 
l(q+ 1)!, and then all (q + 2) columns are scaled by powers 
of /?'//?„ to account for any change in the step size. 

The solution values at the prescribed output times ; ouU , 
f ou , ... are obtained from the Nordsieck history array by 
using the Taylor series expansion method (eq. (A33)) described 
for EPISODE and LSODE. The same procedure used in these 
two codes to start the integration is used in GCKP84; the step 
size, h {) , to be attempted on the first step must be supplied by 
the user. 


CHEMEQ 


In this technique, developed by Young and Boris (ref. 10), 
the species rate equation (Al) is expressed as a difference 
between two positive-definite terms as follows: 


( ~T=f = i= (A3 8) 

dt 

where, for species /, the production rate (P, and the destruction 
rate 3D, can be derived from equation (3): 

N r 

V, = P ' Yj ( y ij' R -j + V "j R :' ] 
j= 1 

Nr 

a = p 1 Y l''j R -j + < R - 

j= i 

When the temperature ODE (eq. (9)) is required (method B), 
it can be cast in a similar form by combining equations (9) 
and (A38) 

N s N s 

- D/a 

Av s +1 _dT _ k=\ _ ^=1 

dr ” dt ~ “ N~s 

D ^ kC p ■ k D ^ kC p ' k 

k = 1 k = I 

= ^ > a , .v+ 1 ~ ^,y v 4- 1> 

= (P 7 — T)y, 



where 


and 


(p, = 


N s 

D ^ 

k= 1 
N k 

D ^ kC P' k 

k — 1 


j Yv 

D 65 A 


D y'k c p.k 
k= 1 

The objective of this decomposition is to enable factorization 
of v, from 3D, 


3D, = £,-v, = V, / r , 

where £ h the loss coefficient for species /, is obtained simply 
by dividing 3D, by v, (i.e., £, = 3D,/v,). With this new 
notation, equation (A38) can be written as 


— = (P, - £;V, = (?i ~ Vi/V (A40) 

dt 


which, for constant (P, and can be solved to give 


(P, (P, 

y.Un) = ViUn- \+K) =— 4 - }'iUn-\) ~ — C*p(- £jh„) 
<L, £j 

Expressed in this way, it can be seen that 1/JC, ( = t,) describes 
how quickly the variable y, reaches its equilibrium value. 

In advancing the solution from time f ;j _ t to time t, n all of 
the equations are separated into two classes, stiff and nonstiff, 
according to the criterion 

h n ( > ] stiff 

j , ( <1 nonstiff 


where r i n _ ] denotes the value of r, at time t n _ x . The two 
types of equations are integrated by separate predictor- 
corrector schemes. For equations classified as nonstiff, the 
improved Euler method (with the Euler method as predictor 
and the modified Euler method (or trapezoidal rule) as 
corrector) is used. For equations classified as stiff, a simple, 
stable, asymptotic formula is used. 
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Nonstiff predictor: 


yw = y ,„_, + /a,,-. 


After each step n the step size h' }+ t to be attempted on the 
next step is computed from the converged integration cycle 
(A41a) as follows: 


Stiff predictor: 

y 10} _ ~ hn) + 2A„ Tj n -)(Pj n - \ (^41b) 


Nonstiff corrector: 


K + 1 = K 


+ 0.005 

(a/EPS ) xu 


The step size, h 0 , to be attempted on the first step is 
determined such that none of the variables will change by more 
than a prescribed amount. The formula used for h {) is 


1 1 = v 

1 i,n 1 i,n — I 



(A42a) 


h 0 = EPS 


mm 


, or — — if V/(^o) 

f(t o) 


< 10 


-20 


Stiff corrector: 


y[m + 1} 
1 i.n 




+ T i,n- 1 




+ 


y , n - 1 



+ 




+ T i.n-\ 



(A42b) 


In equations (A4I) and (A42), m + 1 is the current iteration 
number. The zeroth estimate is the result of the predictor step. 
Also, />< =f(\Yl'"!\)- Convergence is ascertained by com- 
paring Yi”' +]] with Vi” 1 } for all N components using the relative 
error criterion 


max 


I n 


m+ 1 1 _ y ! IW I | 


in Ik, 1 ;:', n'r 11 ! 


min 


9 

< EPS 


(A43) 


To avoid numerical difficulties with the use of equation (A43), 
each estimate is constrained by a minimum value. In the 
present work, a variable that is less than 10 2() is set equal 
to 10 -20 . Thus, for a species with decaying concentration, 
convergence is obtained trivially once Y- t < 10 ~ 20 , and its 
equation is decoupled from the equation set. 

If convergence is not achieved after ITMAX iterations, the 
step length is halved and the step repeated. In this study, a 
value of ITMAX = 5 was used because it minimized the 
execution time for both test problems (ref. 2). If the corrector 
converges after M iterations (M < ITMAX), the step is 
accepted as successful, and the solution is updated 

Y,,, = Y}W /=! N. 


The solution at each output station / {>ut was computed by 
linear interpolation between the computed approximations at 
t n _ | and t n , where t n is the first mesh point that is >/ out : 

y« 0 u«) = k„_, + ,oul ~ (Y„ - y„_,) 

A, * 


CREK1D 

In CREK1D, attention is paid to the distinguishing physical 
and computational characteristics of the induction, heat release, 
and equilibration regimes (refs. 11 to 14). This code consists 
of two algorithms developed for the two distinctly different 
regimes: (a) induction and early heat release, when the ODE’s 
are dominated by positive time constants and (b) late heat 
release and equilibration, when the ODE’s are more stable 
(ref. 2). Both algorithms are based on an exponentially fitted 
trapezoidal rule, but they use different iterative methods for 
convergence. 

The code CREK1D solves a mixed differential-algebraic 
system of equations: ODE’s for the species mole numbers 
and the algebraic enthalpy conservation equation (8) for the 
temperature. The ODE’s and algebraic equation are solved 
simultaneously; however, in the following discussion the 
variables y and Y refer only to the species mole numbers. 

The solution method used for the species ODE’s is a gener- 
alized, tunable, single-step, implicit procedure: 

Y iit1 = Y in _] + h n [UiJ U n + O - V Ln )f,n - iJ / = 1 Afe 

(A44) 

where U i n is a degree-of-implicitness factor. This parameter 
is obtained by “exponentially fitting 1 ' it to a locally exact 
solution of equation (Al) as follows: The species rate expression 
f is expressed in a locally linearized form such that 


No attempt is made either to estimate or control the local 

truncation error. f ~ f.n- 1 — * — 1 -• (A45) 
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where the choice of 0, ,,, a suitable linearization constant, is 
discussed shortly. Equation (A45) assumes that in the interval 
/„] (i.e., locally) each species mole number varies expo- 
nentially. Integration of this equation gives 


iteration in regime b and Jacobi-Newton (JN) iteration (ref. 29) 
in regime a, 

A Newton-Raphson functional F,]”' 1 (/ = 1 N s ) for each 

species mole number is defined from equation (A44) by 


Yi.n - Yi.n- I + K f.n- 


exp (M/.») ~ 1 




(A46) 


Z7 1" 

1 m 


\m\ 


vM _ v 

1 i.n 1 i,n - 

^tAi.n 


i - q , H 

Ui.„ 


An -An" 


(A51) 


To exponentially fit U i n we first replace f n in equation (A44) 
with the expression obtained from equation (A45): 


for i = 1 For temperature the functional Fj n " is 

defined from the enthalpy conservation equation (8) as 


f,n ~f.it— i + 0',n(Y iin - (A47) 

and then equate the resulting expression for Y t n with equation 
(A46). These operations give 


F\"l! = £ Yffl - H 0 (T 0 ) 


(A52) 


k = 1 


U, 


1 1 

_j_ 

h„e, M 1 - exp (h,fi, n ) 


i = 1 N s (A48) 


In order to maintain absolute ^-stability of equation (A48) 
(i.e., to keep errors introduced into the numerical solution at 
any one step bounded as h n is increased indefinitely), U, n 
must be restricted to the interval (0.5, 1.0). For values of 
0, „ > 0, equation (A48) gives U i n <0.5. CREK1D resolves 
this problem by setting 6 i n = 0 whenever it is greater than zero. 
This value of d in gives U i n = 0.5, so that equation (A44) 
defaults to the second-order-accurate trapezoidal rule. However, 
for 8 in < 0, equations (A44) and ( A48) together are equivalent 
to the locally exact or exponential solution, which has an 
equivalent polynomial accuracy of order six to eight (ref. 1 1). 
Thus, equations (A44) and (A48), with the constraint 
(0.5 < U in < 1), constitute an exponentially fitted trapezoidal 
rule, a method which is /I-stable and has a polynomial-order 
accuracy of at least two and as great as six to eight. 

The linearization constants (0 />w j are obtained in one of two 
ways. In the first, called functional linearization (see refs. 1 1 
to 14), equation (A47) is solved explicitly for 0, „ to give 


where m is the iteration number, 7j' f| is the /^-approximation 
to the exact value T{t n ) , it* (7^) is the molar-specific enthalpy 
of species k at temperature T}"’K and H 0 (T 0 ) is the initial 
mixture mass-specific enthalpy at the initial temperature T 0 . 

Newton-Raphson corrector equations with log variable cor- 
rections (for self-scaling of the widely varying mole numbers) 
are given by 


Y Ain >tr" + Ain = -Fft 

U . ain Y[";j ■ din T};"' 


k = I 


for i = 1 N s . 


(A53) 


n k 


[m\ 


my, 

Lj l ain li m J 


Ain Yl‘"„ 


+ 11 


+ 


dF 


\nt j 
r.n 


k= 1 


din T\" a 


Ain 7l”' +l1 = 


" ‘ l\n 


(A54) 


where 


Ain Y};? + " = In K, l '' ,+ I1 


in Yl", 


i = 1 N s 


q fi,n Sm — 1 2 


V — V 

J i,n I i.n— 1 


(A49) 


and 


Ain 7^ +j| = In 7j wl+l1 - In 7j/ w| 


In the second approach, called formal linearization (refs. 1 1 
to 14), the net formation rate of each species is expressed as 
a difference between two positive-definite terms, as described 
in the previous section (see eqs. (A3 8) and (A40)). Comparing 
equations (A47) and (A40) gives 

0i,„= i (A50) 


for this procedure. 

At each integration step, equation (A44) must be solved for 
Yj , r The solution is accomplished by Newton-Raphson (NR) 


The partial derivatives in equation (A53) are obtained from 
equation (A5 1) and are as follows (with the step and iteration 
numbers suppressed for clarity): 


din Y k * \hU, dY k ) 


(A55) 


dF, _ T df, 

din T dT 



here d, k , the Kronecker symbol, is 

= 0 i*k 

= 1 i =k 


The JN iteration technique can be derived from the NR 
iteration procedure by neglecting the off-diagonal elements of 
the Jacobian matrix for the mixed differential-algebraic system 
of equations. With this simplification, equations (A53) to (A55) 
reduce to 


and df/dY k and d//dF can be derived from equations (3) to 
(8), and (10). In evaluating df/dY k , the partial derivatives 
with respect to a m are assumed to be negligible compared 
with the other terms. The required partial derivatives are then 
given by 


an 




Nk 

'E 

7=1 


(»ii 


- v kj R -J 


(A56) 


ain t, 1 ;;" 


Ain Y }"' + 1 1 = -F,!;," 1 i = 1 N s (A59) 


dtr * - Ain 4"' +u = -F, 1 ';; 1 

ain T)l" ] 

yM If |m) 

or i,n _ 1 In _ y |»/1 

ain Y}f~h„u,. n an 1 ;;" 


(A60) 
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dT T 


\ E K _ 

PT “ 


R, 


T 

N, + -*■ 
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N-j + 



The iteration procedure is further simplified by approximating 
dfj/dYj, (eq. (A56)) as follows: 


V, 

an 


(pn> 


n* 

E 

7=1 


( p-j Rj + Vj” R _j) 


where 


which, when combined with equation (A38), gives 



/ = l 



an n 


The partial derivatives of Fj are obtained by differentiating 
equation (A52) and are as follows: 


With these simplifications equation (A59) can be solved 
explicitly for the iterative corrections 


dF r 

ain Y k 




Ain = - 


/rM 

1 iji 

Y):‘}/h„u,,, + 3D);;: 1 


/ = 1 N s ( A6 1 ) 


dFj 

ain T 


N s 


= D 


(A57) 


where, again, the step and iteration numbers have been 
suppressed. The N$ + 1 equations (A53) and (A54) are 
solved simultaneously by LU decomposition and back- 
substitution (ref. 27). The resulting log variable corrections 
are used to update the current estimates (yj” +l1 ! and 7},'" + 1 1 
by the approximate equations 


Y-™ +i] = Y}’:' (l + Ain T,-);; 1 + 1 *) i = 1 N s 

jim+w _ 7i»'l^i + A | n 7j»+ | l'l 


(A58) 


The temperature correction is obtained by substituting equation 
(A57) into equation (A60): 


Ain 7l" ,+ l1 = 


r I\n 


(A62) 


N, 


4"" £ ft;, 1 <-,,*( n m| ) 


* = 1 


where ( 7}/ wI ) is the constant-pressure molar specific heat 
of species k at temperature 7j m| . The current estimates are 
updated by using equation (A58). 

To start this iteration process, the predicted values for the 
species mole numbers are given by equation (A46): 


The solution procedure does not use a predictor; instead, the 
converged results \Y tJl . |j and T n _ { from the previous step are 
used to start the iteration. 


Y] 0 f }=Y Ut .^hJ, n _ 


exp (h n e im „) - 




i = l,..., N s 

(A46a) 
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The predicted temperature is obtained by a single NR iteration 
of the enthalpy conservation equation (8) 


HM - £ YfS Mr,,.,) 

n 01 = t „- i + — - — — 

J \V 

2 ( )i.k(T„- j ) 

A - 1 


(A63) 


If convergence occurs on the first iteration, C m is set equal 
to 0. 1 . 

At each step an average weighted local truncation error 
estimate, d,„ is computed by using the approximations 



1 A,s ( Y — 

1 ^ \ 1 i.n x i.n 

Ns £ 1 (_ max <r,,n.r„) 



For both NR and JN iteration schemes the test for 
convergence of the estimates is based on the magnitudes 
of the log variable corrections, and is given by 


N, 


I :'2 


2 < A|n Y '.» > 2 


/= l 




< EPS (A64) 


This test is used only with variables whose magnitudes are 
greater than 10 ” 20 ; that is, the summation does not include 
species with mole numbers <10~ 2 °. At each iteration the 
estimated convergence rate, C,„, defined as 


is also computed. If convergence is not obtained after ITMAX 
iterations, where ITMAX is the user-supplied maximum number 
of corrector iterations to be attempted, or if C,„ > 0.8, the 
step length is decreased. The new step length is calculated as 
follows: 


for the JN iteration, and 



for the NR iteration. The above summations include only 
species whose mole numbers are greater than 10 2l) . For both 
iteration techniques the step length, /? accv , that would exactly 
satisfy the user-specified local relative error tolerance, EPS, 
is calculated from 


= h H (EPS/d„) ,n 

The step length h,[ +] to be attempted on the next step is 
taken to be 


= min(/? ltcr , h a ,, r 1 0/i„ ) (A65) 

However, if convergence difficulties forced a reduction in the 
step length on the current step, h'+ \ is restricted to 

h „+ 1 — min(/i„, /»,!+ i) (A66) 


/?„ — h n min [0.5, max(0.1, 0.5/C,,,)) 

and the step is retried with the new step size. At least two 
iterations are required to define C,„; on the first iteration <5,,,^ , 
is set equal to 10. A value of ITMAX = 10 was used for both 
problems examined in this study. 

If convergence is achieved in M iterations (M < ITMAX), 
the step is accepted as successful, and the solution is updated: 

x,, = W < = i n s 

T — T l' V l 
1 n 1 n 

After corrector convergence, the step length, /7 i{cr , that 
would produce a convergence rate in the range (0.4,0. 5) is 
estimated as follows: 

A iter = A»(0.4/C m ) l/2 C m < 0.4 

=h„ 0.4 < C„, < 0.5 

=//„(0.5/C m ) 1/2 C„, >0.5 


to prevent a recurrence of the problem. 

CREK1D automatically selects the linearization method and 
the iteration scheme to be used for solving equation (A44). 
During induction and heat release, when small step lengths 
are required for solution accuracy, the JN iteration is used 
to minimize computational work. During late heat release and 
equilibration, when the ODEN are more stable and larger step 
lengths can be used, the NR iteration is preferred since it has 
better convergence properties than the JN iteration. The regime 
identification test exploits the fact that during equilibration 
many reactions achieve a condition in which the forward and 
reverse rates are large but with vanishingly small differences 
(refs. 13 and 30). The actual test used at the beginning of each 
time step is 

9 

\fi\ < 10 ((P, + 3D,) (A67) 

where (P, and 2D, are given by equation (A39). If any two 
species satisfy equation (A67), regime b is obtained, and the 
NR iteration is used for the step. If fewer than two species 
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satisfy equation (A67), regime a is obtained, and the JN 
iteration is used for the step. Once the NR iteration is selected 
for any one step, the above test is no longer applied, and the 
NR iteration is used for the rest of the problem. 

Whenever the reaction rate for any species satisfies equation 
(A67), that species is considered to be in “quasi-steady state" 
and the “L-formulated" equation (A50) is used. For all other 
species the “Z- formulated" equation (A49) is used. To minimize 
computational work, the [Z i n j are evaluated only once per 
step: at the beginning of the time step, using equation (A49). 
However, since Y { „ and/ ,, are not known at the start of the 
step, the \Z jn j are approximated using values from the 
previous step: 

„ y fi,n — 1 /,» — 2 

A.h — A.w- I — v _ v 

M./i— 1 */,« — 2 

CREK1D also includes an algorithm for filtering the initial 
conditions that may be ill posed. These ill-posed conditions 
may arise, for example, in multidimensional modeling because 
of the averaging of mole numbers over adjacent grid nodes. 
CREK ID therefore “filters" the initial conditions to provide 
physically meaningful initial mole numbers and net species 
production rates. For purposes of this filtering CREK ID uses 
the decomposition performed in CHEMEQ (eqs. (A38) and 
(A40)). On the first call to CREK ID it uses this formulation 
over one time step of length A], given by 

A, = (A68) 

max <£/ o 

i 

The predictor-corrector algorithm uses equation (A46) (with 
6, , = -<£, o) as the predictor 


Equation (A69) is iterated until converged; that is, the criterion 
given by equation (A64) is satisfied. It convergence is not 
obtained after 10 iterations or C„, > 0.8, the step length is 
halved, and the step retried. If convergence is obtained after 
M iterations (M < 10), the step is accepted as successful, the 
solution for the mole numbers is updated 

Y U = Y}P i=l Vs- 

and the temperature T\ is obtained by a single Newton- 
Raphson iteration 

Av 

HM - Yt Y k.MT») 


£ Yk. .<>.*( r 0 > 

A = 1 

The step size, to be attempted on the next step is 
determined from the maximum loss coefficient at t\ by using 
an expression similar to equation (A68). For this step, the JN 
iteration (eqs. (A46a) and (A61) to (A63)) is used, with all 
9; n set equal to zero, so that all U x n = 0.5 (see eq. (A48)). 
The predictor step (eq. (A46a)), therefore, reduces to the 
explicit Euler method, and the corrector (eq. (A44)), to the 
trapezoidal rule. For the next and subsequent steps the step 
size is adjusted according to equation (A65) or (A66), and the 
iteration procedure and linearization constants are selected as 
described previously. If NR iteration is used, the Jacobian 
matrix for the mixed differential-algebraic system of equations 
is updated at y = Y n ~\, T—T n _\. 

The solution values at the prescribed output times ? out j, 
r oul2 ,... are obtained by adjusting the step length so that the 
internal mesh points coincide with these times. Thus, the step 
size A,; + | is given by 




i= 


An implicit Euler corrector is then iterated to convergence 

tir" = r i .o + hdr + " 

In the above two equations, the subscript 1 indicates that this 
is the first step. Using equations (A38), (A40), and (A58), 
together with the approximations (Pj [ +!} = (P}"} 1 and + 1 1 = 
£)'i , l, the preceding corrector equation can be rewritten to 
provide the following expression for the log variable 
corrections [Ain Y$ +l *\: 


Ain Yl;i ,+ l] = 


Yj . o - Yl T 1 + A ,/ 1 ? 1 
Yl r 1 + A,3Djj 3 


i=\,...,N s (A 69) 


K+\ - min r oul -t„), 

where f oul is the current value of the output time, and the 
results at r out are generated by solving the governing 
equations. To continue the integration past each output time, 
the procedures described above for the second and subsequent 
steps are used. 

To reduce the computational cost, the use of exponential 
functions is minimized by replacing them w ith rational function 
approximations. For example, the term (e x — l)/ar in equation 
(A46a) is evaluated by means of a (2,2) diagonal Pade 
approximation, e x (2 , 2 )^ f° r ex P x: 


e b,2) 


x 

1 + - 
2 



x x~ 
— + — 

2 12 


x < 0 
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which gives 

1 < 0 

, 1 - v(l/2 -.v/12) 

Similarly, the tuning factor U i<n (cq. (A48)) is evaluated by 
using the approximation 



This equation requires six operations to evaluate and does not 
exhibit the singularity at x = 0 of the exact expression (eq. 
(A48)). 


Although log-variable corrections are used in the code, 
evaluation of logarithms of the variables is not required. Also, 
the use of the approximations given by equations (A58) avoids 
the cost of computing the exponentials of the log-variable 
corrections to obtain the new estimates. 

Another technique used in the code to reduce computational 
work is to locally linearize the expressions for the thermo- 
dynamic properties of the species and the rate coefficients. 
In particular, during the course of iterative convergence of 
the equations, the thermodynamic properties and rate 
coefficients are not reevaluated while the current temperature 
is within a local window (T, T + AT ) , where AT is specified 
by the user. Use of this strategy has been shown to reduce 
the computational work (refs. 2, 3, and 5). 
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Appendix B 

Description of Temperature Calculation Method A 


In this method, the temperature T„ at each discrete time 
is computed from the solution for the species mole numbers 
returned by the integrator by using the algebraic enthalpy 
conservation equation 


As 

£ o,',MT n ) = H n (8) 

i= 1 

Equation (8) is solved for the temperature by using the 
Newton-Raphson (NR) iteration technique (e.g., ref. 23). This 
equation is rewritten as 


Ns 

F(T„)= £ a,J,(T„) - H 0 (Bl) 


is given by 



< ERMAX 


where the vertical bars denote absolute value and ERMAX 
is the local relative error tolerance. If convergence is not 
obtained after MAXITS iterations, where MAXITS is the user- 
supplied maximum number of corrector iterations to be 
attempted, an error exit is taken. If convergence is achieved 
in M iterations (M < MAXITS), the solution T}, m * is accepted 
as the temperature at 


T — j\ M\ 
1 n 1 n 


so that solving equation (8) is equivalent to finding the zero 
of F. The quantity F(T}”^) is the amount by which the mixture 
mass-specific enthalpy at the m ih approximation for T n , T}” ,] 
(m = 1,2,...), fails to satisfy equation (8). A new approximation, 
T}" l+{ K for the temperature is obtained from equation (Bl) 
by locally linearizing F at T} t m h 


7*»i+i| _ 7i»i| __ ^ m j 

" " (dF/dT) T=T m 


i= 1 


The test for convergence of the estimates is based on the mag- 
nitude of the corrections 67j, w+l * (equal to 7]| /N+l1 — T}^) and 


The NR iteration will converge if the initial guess (i.e., 
n 01 ). is sufficiently accurate (ref. 23). The present work did 
not utilize a predictor to generate 7j 01 ; instead, the most recently 
computed temperature was used to start the iteration. Now, 
the temperature was evaluated at the end of each integration 
step and whenever the species derivatives and Jacobian matrix 
were computed. Hence, the converged value obtained either 
at the end of the previous step or from the previous estimates 
for the mole numbers was used as the initial guess for the current 
temperature. For the very first temperature computation for 
the problem, the initial temperature, 7J), served as the predicted 
value. The above procedure was found to be satisfactory in 
that the iteration converged for all integration methods and 
EPS values used in this study. In addition, the converged temper- 
ature was not significantly different from that obtained by 
integrating the temperature differential equation (ref. 2). 
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